Evolution of warped and twisted accretion discs in close binary systems
aa r X i v : . [ a s t r o - ph . E P ] D ec Astronomy&Astrophysicsmanuscript no. manusscript13088 c (cid:13)
ESO 2018November 8, 2018
Evolution of warped and twisted accretion discs in close binarysystems
Moritz M. Fragner & Richard P. Nelson.
Astronomy Unit, Queen Mary, University of London, Mile End Road, London E1 4NS.e-mail:
[email protected], [email protected]
Received / Accepted
ABSTRACT
Context.
There are numerous examples of accretion discs in binary systems where the disc midplane is believed to be inclined relative to thebinary orbit plane.
Aims.
We aim to examine the detailed disc structure that arises in a misaligned binary system as a function of the disc aspect ratio h , viscosityparameter α , disc outer radius R , and binary inclination angle γ F . We also aim to examine the conditions that lead to an inclined disc beingdisrupted by strong di ff erential precession. Methods.
We use a grid-based hydrodynamic code to perform 3D simulations. This code has a relatively low numerical viscosity comparedwith the SPH schemes that have been used previously to study inclined discs. This allows the influence of viscosity on the disc evolution to betightly controlled.
Results.
We find that for thick discs ( h = .
05) with low α , e ffi cient warp communication in the discs allows them to precess as rigid bodieswith very little warping or twisting. Such discs are observed to align with the binary orbit plane on the viscous evolution time. Thinner discswith higher viscosity, in which warp communication is less e ffi cient, develop significant twists before achieving a state of rigid-body precession.Under the most extreme conditions we consider ( h = . α = × − and α = . ff erential precession. Discs that become highly twisted are observed to align with the binary orbit plane on timescales much shorterthan the viscous timescale, possibly on the precession time. Conclusions.
We find agreement with previous studies that show that thick discs with low viscosity experience mild warping and precessrigidly. We also find that as h is decreased substantially, discs may be disrupted by strong di ff erential precession, but for disc thicknesses thatare significantly less ( h = .
01) than those found in previous studies ( h = .
1. Introduction
Accretion discs are found in a number of astrophysical environ-ments. These include protostellar discs around T Tauri stars,and around compact objects in cataclysmic variable and X-ray binary systems. Most T Tauri stars associated with low-mass star-forming regions are observed to be members of bi-nary systems, where the peak in the distribution of separa-tions occurs at approximately 30 AU (Leinert et al. 1993).Studies of the tidal interaction between a circumstellar discand a binary companion indicate that tidal truncation oc-curs at a ratio of binary separation to disc radius D / R ≃ ff ects.There are observational indications that the disc and binaryorbital plane are not always coplanar. The binary system HKTau has been imaged, and shows compelling evidence for adisc whose midplane is misaligned with the orbit of the bi- nary companion (Stapelfeldt et al. 1998). Terquem & Bertout(1993, 1996) have suggested that the reprocessing of radiationfrom the central star by a tilted and warped accretion disc couldaccount for the high spectral index of some T Tauri stars. Anumber of precessing jets have been observed in star formingregions, and the jet precession may indicate an underlying pre-cessing disc from which the jet is launched (K¨onigl & Ruden1993; Terquem et al. 1999). In the case of close X-ray binaries,the objects SS433 and Her X-1 show evidence of precessingjets and / or accretion discs (Boynton et al. 1980; Katz 1980;Margon 1984; Petterson 1977).An early study of the equations describing the behaviour ofa thin non-coplanar viscous disc was undertaken by Petterson(1977). The resulting equations, however, lacked the propertyof conserving global angular momentum. This was addressedin subsequent work by Papaloizou & Pringle (1983), wherethey showed that warps in discs for which h < α < ff usively on a timescale shorter than the usual viscous ac-cretion by a factor α , where α is the dimensionless viscositycoe ffi cient (Shakura & Sunyaev 1973), and h is the disc aspect M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems ratio H / r . The opposite regime α ≤ h has been investigatedby Papaloizou & Lin (1995) using linear perturbation theory.They demonstrated that the governing equation for the disc tiltchanges from being a di ff usion equation to a wave equation inthis physical regime, showing that warps propagate as bendingwaves with a speed corresponding to half the sound speed, c s / ff erential precession time. Furthermore, they estimate thatthe timescale for damping the inclination is of the same order asthe accretion timescale of the disc. A study of the twisting andwarping of viscous discs in binary systems was undertaken byLubow & Ogilvie (2000) using linear theory, and it was shownthat discs whose outer disc radii are smaller than the truncationradius maybe unstable to tilting.Numerical simulations of tidally interacting discs whichare not coplanar with the binary orbit have been performed byLarwood et al. (1996) using SPH simulations. They found thatthe disc precesses approximately as a rigid body as long as thedisc aspect ratio is large enough. For smaller aspect ratios theirresults suggest that the disc develops a modest warp, while for h ≤ .
03 the disc becomes disrupted as a consequence of dif-ferential precession. Their results suggest that a disc can split-up into two distinct parts that behave independently. TheseSPH simulations also showed that the inclination evolves onthe viscous timescale, as expected from Papaloizou & Terquem(1995).The goal of the present study is to examine in detail thestructure of misaligned accretion discs in close binary systemsas a function of the important physical parameters: h , α , theouter disc radius R , and the inclination angle γ F . The binarycompanion and the central star are assumed to be of equal mass.We use a grid-based code (NIRVANA) to perform this study,and this code has a relatively small numerical viscosity com-pared with the SPH schemes used to perform previous studies.This allows us to tightly control the disc viscosity, and thereforeexamine in more detail its influence on the disc evolution. Weconsider a broad range of disc parameters in which warps prop-agate either as bending waves h > α or di ff usively h < α . Aswell as examining the quasi-steady disc structures which arisebecause of tidal interaction with the inclined companion, wealso examine the conditions under which a disc becomes dis-rupted due to strong di ff erential precession. In particular we areinterested in examining whether the following criterion holds: adisc will achieve a state of rigid-body precession if warp prop-agation across the disc (either through bending waves or di ff u-sion) occurs on a timescale shorter than the di ff erential preces-sion time. The plan of the paper is as follows. In Sect. 2 we presentthe basic equations of the problem. In Sect. 3 we describe thenumerical methods used in the simulations, and in Sect. 4 wediscuss the linear theory of bending waves and calibrate thesimulation code against calculations based on linear theory. InSect. 5 we present our results for the tilted, tidally interactingdiscs, and in Sect. 6 we discuss our results and present our con-clusions.
2. Basic equations
We consider the evolution of a thin, gaseous, viscous disc inorbit around a central star, where the disc is perturbed by acompanion star orbiting in a plane which is not coincident withthe disc midplane. For a range of physical parameters, it is ex-pected that the disc will precess rigidly around the angular mo-mentum vector of the binary system. The equations of conti-nuity and motion for a viscous fluid in a precessing referenceframe may be written as: D ρ Dt + ρ ∇ . v = D v Dt = − ρ ∇ P − ∇ Ψ + S visc − Ω × v − Ω × ( Ω × r ) (2)where DDt = ∂∂ t + v · ∇ is the convective derivative, ρ is the density, v the velocity and Pthe pressure. The precession frequency is given by | Ω | , and thedisc angular momentum vector is assumed to precess around avector pointing in the direction of Ω . The gravitational poten-tial is Ψ . We adopt a locally isothermal equation of state: P = c s ( r ) ρ (3)where the isothermal sound speed is defined by c s ( r ) = h · v Kep .Here h is the aspect ratio of the disc, H / r , and v Kep is the localkeplerian velocity. S visc is the viscous force per unit mass: S visc = ρ ∇ T (4)where T is the viscosity stress tensor: T i j = ρν ∂ v i ∂ x j + ∂ v j ∂ x i − δ i j ∇ . v ! (5)We adopt the standard ‘alpha’ model (Shakura & Sunyaev1973) to specify the kinematic viscosity ν = α c s H . We solvethe above equations numerically, using a system of sphericalcoordinates r = ( r , θ , φ ). We work in a reference frame in which the origin of the coordi-nate system is fixed on the central star, and the secondary star .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 3 moves on a circular orbit about the origin with position vec-tor D . The gravitational potential, Ψ , at any position vector r istherefore given by: Ψ = − GM P r − GM S | r − D | + GM S r · D D , (6)where G is the gravitational constant, and M P and M S are theprimary and secondary masses, respectively. The first term isdue to the primary star, and the remaining two terms are dueto the companion star. The last indirect term accounts for theacceleration of the origin of the coordinate system, which co-incides with the location of the primary star. Note that there areno contributions to the potential due to the disc since we do notinclude disc self-gravity in our calculations, and the stars do notexperience a gravitational force due to the disc. Consequentlythe binary orbital parameters remain unchanged.Hence the secondary star feels the acceleration due to thecentral star and the indirect component of the potential. In anon-precessing reference frame, centred on the central star,whose x − y plane coincides with the orbital plane of the bi-nary, the position vector of the companion in a circular orbitmay be written as: D = D · (cos ( ω t ) e + sin ( ω t ) e ) , (7)where D is the constant binary separation and ω = p G ( M P + M S ) / D is the binary orbital angular velocity.
3. Numerical methods
The system of equations described in Sect. 2 is inte-grated using the grid-based hydrodynamics code NIRVANA(Ziegler & Yorke 1997), adapted to solve the equations in a pre-cessing reference frame. This code uses operator-splitting, andthe advection routine uses a second-order accurate monotonictransport algorithm (van Leer 1977).
Periodic boundary conditions were applied in the azimuthaldirection. At all other boundaries outflow conditions wereadopted using zero-gradient extrapolation. Velocities at the ra-dial boundaries were set to ensure that zero viscous stress oc-curs there.
We adopt a system of units such that the unit of mass is equalto that of the central star, M P , the unit of radius is equal tothe radius of the inner boundary of our computational domain,and we set the gravitational constant G =
1. The unit of timethen becomes 1 / Ω K (1), where Ω K is the keplerian angular ve-locity evaluated at r =
1. When discussing simulation results,however, we will express time in units of the orbital periodat r =
10, the nominal outer disc radius in the simulations.Thus we describe a time interval of 2 π Ω K (10) − as ‘one orbit’.Inclination and precession angles are described in units of de-grees. The numerical domain extends radially from r = [1 , R ], andazimuthally from φ = [0 , ◦ ]. The outer radius, R , may di ff erfrom simulation to simulation. The meridional domain coversthe range θ = [90 ◦ − ∆ θ, ◦ + ∆ θ ], where again ∆ θ may varybetween the simulations.We perform simulations in two di ff erent reference frames.For models in which the disc is expected to develop rigid pre-cession, we solve the equations in a frame which precessesaround the angular momentum vector of the binary. Given asensible choice of the precession rate, this approach ensuresthat the disc midplane always stays close to the equatorial planeof the computational domain where θ = ◦ , and allows sim-ulations to be performed with large inclinations of the binaryorbit. If evolved in a non-precessing frame, the disc would pre-cess around the binary angular momentum vector, causing itto eventually intersect with the upper and lower meridionalboundaries of the computational domain if the binary inclina-tion angle is large. Adopting a precessing reference frame al-lows simulations to be conducted with relatively small merid-ional domains, even for large binary inclinations, thus substan-tially reducing the computational expense involved.For models in which substantial di ff erential precession ofthe disc is expected to arise, adopting a frame with a single pre-cession frequency does not solve the problem described above,since some parts of the disc inevitably intersect the meridionalboundaries. As these models are of significant interest fromthe astrophysical point of view (they will tend to be thin discswith large viscosity, similar to those which occur in X-ray bina-ries for example), we have computed them in a non precessingframe, but have been forced to use large meridional domainsand smaller binary inclination angles.We denote coordinates in the precessing frame as ˆr , whilecoordinates in the non-precessing binary frame are denoted r .The coordinates in the code frame are denoted with r C , wherethe code frame can be one of either the precessing frame orthe binary frame. The transformation of any vector e from thebinary into the precessing frame is given by: ˆe = R Z ( Ω F t ) R X ( γ F ) e (8)with rotation matrices: R X ( γ F ) = γ F ) sin( γ F )0 − sin( γ F ) cos( γ F ) R Z ( Ω F t ) = cos( Ω F t ) − sin( Ω F t ) 0sin( Ω F t ) cos( Ω F t ) 00 0 1 . (9)Thus to transform a vector from the binary into the precessingframe, we first rotate around the x -axis by an inclination an-gle γ F and then rotate the resulting vector around the z -axis byan precession angle Ω F t . Here Ω F is the precession rate of theprecessing frame and t is the time elapsed. This transformationis particularly useful when setting up the initial conditions forthe simulations which are performed in the non-precessing bi-nary frame, since we assume that the binary orbit lies in a planewhich is coincident with the equatorial plane of the computa-tional grid, with the disc being inclined by an angle γ F . M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems
In the precessing frame the code midplane (defined to residewhere θ C = ◦ ) coincides with the disc midplane initially, sothat ( r C ≡ ˆr ). The initial density profile is chosen to satisfy hy-drostatic equilibrium in the ˆ θ -direction, which yields the result: ρ ( r , ˆ θ ) = r − ζ h sin (ˆ θ ) i (cid:16) h − ( ζ + (cid:17) , (10)where ζ = is the radial power law exponent of the densityprofile. Note that this expression has the limit:lim ˆ θ → ρ ( r , ˆ θ ) = r − ζ exp − (cid:16) ˆ θ − π (cid:17) h , (11)which is the familiar form of the density profile for thin discs.The velocity in the ˆ φ direction is chosen to balance the forcesin the radial direction:ˆ v φ = r (1 − ( ζ + h ) . (12)Consequently, the azimuthal velocity profile is slightly subke-plerian because of the radial pressure gradient. The velocitiesin the r and ˆ θ direction are chosen to be zero initially.We set Ω in Eq. (2) to point in the same direction as the binaryangular momentum vector: Ω = Ω F e = Ω F ( − sin ( γ F ) ˆe + cos ( γ F ) ˆe ) (13)where the last equality is derived using the transformationgiven by Eq. (8). The precession rate may be estimated usingthe expression: Ω F = − GM S D ! R R Σ r dr R R Σ r Ω dr cos ( γ F ) , (14)where Σ = R ∞−∞ ρ r d θ is the disc surface density. Hence the pre-cession frequency is determined by the ratio of the integratedexternal torque exerted by the secondary to the integrated an-gular momentum content of the disc (Papaloizou & Terquem1995). For the velocity and density profiles given by Eqs. (12)and (10), respectively, this evaluates to: Ω F Ω d ( R ) = − R D M S M P cos ( γ F ) , (15)where Ω d ( R ) is the orbital angular velocity of the disc at radius R . The binary orbit given by Eq. (7), expressed in precessingframe coordinates, is given by: D = D cos ([ ω − Ω F ] t )sin ([ ω − Ω F ] t ) cos( γ F )sin ([ ω − Ω F ] t )) sin ( γ F ) . (16)Thus an observer moving with the disc sees an increased binaryfrequency ω − Ω F due to the retrograde precession of the disc(i.e. Ω F is negative). When working in the non-precessing binary frame we set theprecession frequency Ω = r C ≡ r ), the disc is inclined with respect to theequatorial plane of the computational grid by an inclination an-gle γ F . We can use the transformation defined by Eq. (8) at t = (ˆ θ ) = cos ( φ ) sin ( θ ) + cos ( γ F ) sin ( φ ) sin ( θ ) + sin ( γ F ) cos ( θ ) − ( γ F ) sin ( γ F ) sin ( θ ) cos ( θ ) sin ( φ ) (17)which can be used to determine the initial density profile givenby Eq. (10). The initial velocity is given by: ˆv φ = ˆ v φ ( − sin ( ˆ φ ) ˆe + cos ( ˆ φ ) ˆe ) = ˆ v φ ( − sin ( ˆ φ ) e + cos ( ˆ φ ) cos ( γ F ) e − cos ( ˆ φ ) sin ( γ F ) e )Using the relationship between the polar angles one can findan expression for ˆv φ that is written entirely in terms of binaryframe angles. The di ff erent components in the radial, azimuthaland meridional direction are then given by: v r = v φ = ˆ v φ cos ( γ F ) sin ( θ ) − sin ( φ ) sin ( γ F ) cos ( θ )sin (ˆ θ ) v θ = ˆ v φ sin ( γ F ) cos ( φ )sin (ˆ θ ) (18)where ˆ v φ is given by Eq. (12) and sin (ˆ θ ) by Eq. (17). It canbe shown that this set of profiles satisfies the Euler equationsgiven by Eqs. (1) and (2) with S visc = M S = Ω = γ F = constant, since in a spherically symmetric potential thereis no preferred direction for the disc rotation axis. In order to follow the evolution of the disc structure we calcu-late the total angular momentum vector, L ( r ), for disc materialcontained in each spherical shell in the numerical domain: L ( r ) = Z π d φ Z π +∆ θ π − ∆ θ l · r sin ( θ ) d θ δ r (19)where δ r is the radial grid spacing and l is the angular momen-tum density: l = ρ r × v = ρ ( − rv φ θ C + rv θ φ C ) . (20) θ C and φ C are the code unit vectors in the meridional and az-imuthal directions, respectively. Note that when working in theprecessing frame a term involving the precessional velocitiesshould be added to this expression. In practice we find that thisis of negligible magnitude, because of the slow rate of preces-sion, and therefore the term has been omitted. For disc material .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 5 located in a given spherical shell we locate the inclination an-gle, δ , (equal to γ F at time t =
0) between the disc and binaryangular momentum vectors through:cos ( δ ) = J B . L | J B || L | . (21)where J B represents the binary angular momentum vector. Theprecession angle, β , measured in the binary plane can be de-fined through:cos ( β ) = ( J B × L ) . u | J B × L | , (22)where u is an arbitrary reference unit vector in the binary plane.We choose u = − e , so that the initial precession angle β = − e hasto be expressed in precessing frame coordinates and becomestime dependent.When discussing the results of the simulations we use thethe following terminology: a warp occurs if the angle δ be-comes a function of radius in the disc; a disc develops a twist when β becomes a function of radius; a disc is said to be bro-ken if either β or δ change discontinuously at some radius suchthat the disc separates into two independently precessing parts;a disc is said to be disrupted if β varies smoothly by more than180 degrees across the disc radius because of di ff erential pre-cession. Later in this paper we will be interested in measuring the dif-ferent contributions to the changes in the precession and incli-nation angles for disc material located at di ff erent radii in thecomputational domain. As a first step we calculate the rate ofchange of the angular momentum vector L ( r ) associated withdisc material located at radius r , which involves integratingover the individual spherical shells in the computational do-main.Using the continuity and momentum Eqs. (1) and (2) toexpress the velocity and density changes in Eq. (20), and as-suming that the density vanishes at the meridional boundaries,we find the result: ˙L ( r ) = ˙L F ( r ) + ˙L V ( r ) + ˙L S ( r ) (23)where a dot denotes a time derivative. We have that ˙L F ( r ) = − Z π d φ Z π +∆ θ π − ∆ θ d θ sin ( θ ) h r · l · v r i r + δ r r − δ r (24)is the change due to the divergence of a radial angular momen-tum flux. The change due to viscous interaction with neighbor-ing shells is given by: ˙L V ( r ) = Z π d φ Z π +∆ θ π − ∆ θ r sin ( θ ) d θ h − T r φ θ C + T r θ φ C ) i r + δ r r − δ r . (25)Thus the z -component of the angular momentum vector canonly be changed by the T r φ component of the viscous stress ten-sor, while the x and y components can also be changed due toradial shear of vertical motion via T r θ . The last term in Eq. (23) which induces a change of the angular momentum vector ofdisc material contained in a spherical shell is due to the gravi-tational interaction with the secondary star: ˙L S ( r ) = ˙L T ( r ) + ˙L A ( r ) , (26)where the angular momentum change due to torques is givenby: ˙L T ( r ) = Z π d φ Z π +∆ θ π − ∆ θ t · r sin ( θ ) d θ δ r (27)and the torque density t is given by: t = − ρ r × ∇ Ψ = ρ " θ ) ∂ Ψ ∂φ θ C − ∂ Ψ ∂θ φ C . (28)When working in the precessing frame there is an additionalterm causing angular momentum change due to Coriolis andcentrifugal forces: ˙L A ( r ) = Z π d φ Z π +∆ θ π − ∆ θ ρ r sin ( θ ) d θ (cid:16) − r ∆ v φ θ C + r ∆ v θ φ C (cid:17) δ r with ∆ v φ = ( − Ω × v − Ω × ( Ω × r )) φ C ∆ v θ = ( − Ω × v − Ω × ( Ω × r )) θ C . (29)In the following we are only interested in the torque compo-nents expressed in precessing frame coordinates. Thus all thetorque components have been calculated now for simulationsperformed in the precessing frame. However, when working inthe binary frame, we perform the above calculations and thenproject each of the torque components onto precessing framecoordinate axes to ensure a uniform approach to monitoringthe torque contributions. Since the latter are time dependent,the additional term when working in the binary frame is givenby: ˙L A ( r ) = − Ω × ˆL This term accounts for the angular momentum change mea-sured in the precessing frame that arise because of the preces-sion of the frame.
We can now relate the above torques to the rate of change of theprecession and inclination angles. When working in the pre-cessing frame, the inclination angle δ and precession angle β for disc material in a given radial shell, as defined by Eqs. (21)and (22), are given by:cos( β ) = cos( Ω F t )(sin( γ F ) L Z + cos( γ F ) L Y ) + sin( Ω F t ) L X h (sin( γ F ) L Z + cos( γ F ) L Y ) + L X i cos( δ ) = − sin( γ F ) L Y + cos( γ F ) L Z | L | . (30)Note that the vector components L X , L Y and L Z in the aboveexpression, and in the remaining expressions in this subsec-tion, should have hat symbols because they are measured in M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems the precessing frame. We neglect these, however, in order tosimplify the notation. For su ffi ciently rigidly precessing discsthe deviation from the mean precession frequency Ω F will besmall, and the disc midplane will remain close to the equatorialplane of the computational domain. In this case we can use theapproximation: L X ≪ L Z L Y ≪ L Z . (31)With this approximation, Eq. (30) gives to first order: δ = γ F + L Y L Z β = Ω F t − γ F ) L X L Z . (32)Thus for su ffi ciently rigidly precessing discs the total inclina-tion is the sum of the inclination of the precessing frame and asmall deviation term which is proportional to the y -componentof the angular momentum vector, measured in the precessingframe. Likewise the precession angle is the sum of precessiondue to the frame and a small deviation term that is proportionalto the x -component of the angular momentum vector. The ad-vantage of this approximation is that the precession and inclina-tion angles become linear functions of the torque components.This allows us to write: ∂δ∂ t = ∂δ∂ t ! F + ∂δ∂ t ! V + ∂δ∂ t ! S ∂β∂ t = Ω F + ∂β∂ t ! F + ∂β∂ t ! V + ∂β∂ t ! S (33)where each of the rates are calculated according to Eq. (32): ∂δ∂ t ! K = L Z ˙ L KY − L Y L Z ˙ L KZ ∂β∂ t ! K = − γ F ) L Z ˙ L KX − L X L Z ˙ L KZ . (34)Hence, in circumstances where the approximation expressedin Eq. (31) holds, we can measure separate contributions tothe inclination and precession rate coming from the radial flux(Eq. 24), viscous stress (Eq. 25) and gravitational interactionwith the secondary star (Eq. 26). When discussing contribu-tions to di ff erential precession rates later in the paper, we willomit the constant mean precession rate Ω F of the precessingframe in the expression for the total rate given by Eq. (33).
4. Linear theory of free warps and bending waves
In order to calibrate the code, we now examine how wellit is able to model the propagation of free bending waves,and compare the results with linear theory. The linear the-ory of warps in accretion discs was initially investigatedby Papaloizou & Pringle (1983). Their analysis is valid ina regime in which the dimensionless viscosity parameter α (Shakura & Sunyaev 1973) lies in the range h ≤ α ≤
1. Inthis regime, globally warped discs were found to evolve dif-fusively, with a di ff usion coe ffi cient larger than that associ-ated with mass flow through the disc by a factor of ∼ α , assuming an isotropic viscosity. Papaloizou & Lin (1995) in-vestigated the regime in which α ≤ h , and showed under that,for low frequency modes, the governing equation for the disctilt changes from being a di ff usion equation to a wave equation,such that warps are communicated through the disc via propa-gation of bending waves. In a keplerian disc these are expectedto propagate with little dispersion at half the sound speed, c s / ffi cult to control and spec-ify ab initio .In this section we compare our numerical results witha linear calculation, adopting a similar set-up to the oneused by Nelson & Papaloizou (1999). The equations describ-ing the evolution of linear bending waves have been derivedby Nelson & Papaloizou (1999). We do not reproduce the fullderivation here, but rather mention some of the salient pointsassociated with the derivation. Small-amplitude warps can be considered to be linear pertur-bations of a disc with midplane initially coincident with the( e , e ) plane. Taking the φ dependence of the perturbationsto be ∝ e im φ , with m = r , φ, z ) can be linearised(Papaloizou & Lin 1995; Nelson & Papaloizou 1999). An ana-lytical solution to the linearised equations corresponds to thecase of a rigid tilt. In this case the velocity and pressure pertur-bations are given by: v ′ r = − z Ω K δ v ′ φ = − iz δ d ( r Ω K ) drv ′ z = r Ω K δ W ′ = ρ ′ c S ρ = − irz Ω K δ (35)where Ω K is the Keplerian angular frequency and δ is a smallconstant inclination. For large scale warps, it is expected thatthe local inclination δ varies on a length scale significantlylarger than the disc thickness H = h · r . It is then reasonableto assume the inclination is also approximately independent of z . Using these assumptions, and the rigid tilt solution (35), the z dependence of the linearised equations disappears when in-troducing the new quantities q ′ r = v ′ r z and q ′ φ = v ′ φ z . The resultingexpressions are given by Nelson & Papaloizou (1999) in theform: ∂ q ′ r ∂ t + i Ω K q ′ r − Ω K q ′ φ = i ∂ ( r Ω K δ ) ∂ r ∂ q ′ φ ∂ t + i Ω K q ′ φ + q ′ r r ∂ ( r Ω K ) ∂ r = − Ω K δ .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 7 Fig. 1.
Time evolution of free warp for h = δ MAX = ◦ . NIRVANA results (blue solid line); linear results (red dashed line). ∂ v ′ z ∂ t + i Ω K v ′ z = ir Ω K δ F r Ω K ∂δ∂ t + i Ω K δ ! = − ir ∂ ( r µ q ′ r ) ∂ r + µ q ′ φ r + i Σ v ′ z (36)where F = Z ∞−∞ ρ z c S dz µ = Z ∞−∞ ρ z dz Equations (36) are a set of coupled di ff erential equations whichdescribe the evolution of the inclination δ as a function of r and t . Thus they can be used to follow the evolution of a bendingdisturbance from initial data. Below we use the time-dependentsolution calculated in this way to test NIRVANA’s ability topropagate bending disturbances.Assuming a time dependence of the perturbations ∝ e i σ t ,with the wave frequency obeying the relation σ ≪ Ω K forslowly varying warps, one can derive algebraic equations forthe perturbed quantities. These show that warps induce hor-izontal motions and vertical shear in the disc through termsproportional to (cid:16) ∂ v ′ r ∂ z , ∂ v ′ φ ∂ z (cid:19) (Papaloizou & Lin 1995). Physically,these horizontal motions are generated by pressure forces in thedisc which arise from the vertical misalignment of neighboringregions of the disc midplane. They are responsible for drivingthe bending wave, which propagates with approximately halfthe sound speed.The main e ff ect of viscosity in the disc is to damp these hor-izontal motions, reducing the e ffi cacy of bending wave prop- agation. As the viscosity increases and α > h , the bend-ing disturbance begins to evolve di ff usively on a time scalethat scales inversely with ν/ (2 α ) (Papaloizou & Pringle 1983),where this latter quantity is now the di ff usion coe ffi cient asso-ciated with warp propagation in this regime. If the viscosity issmall enough that its e ff ect on the unperturbed quantities canbe neglected, then the set of equations (36) can be extendedto include the e ff ect of a small viscosity by the replacement(Papaloizou & Lin 1995): ∂∂ t + i Ω K ! → ∂∂ t + ( i + α ) Ω K ! . (37)If the viscosity is large, however, then the complete viscousstress tensor should be included in the equations of motion,making the analysis rather complicated, since there would thenbe accretion velocities in the unperturbed state. A number of simulations using NIRVANA have been per-formed, in which warps with di ff erent amplitudes in discs withdi ff erent parameters were set up and allowed to evolve. Wecompare these solutions with those obtained using the lineartheory expressed in Eqs. (36) using the same set of initial con-ditions. The linear calculations employed a one dimensional fi-nite di ff erence scheme for which we used 1000 equally spacedgrid cells. Starting with a gaussian initial inclination profile: δ ( r ) = δ MAX e ( −| r − r MID | ) M.M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems
Fig. 2.
Final stage of evolution for di ff erent sets of disc parameters. NIRVANA solution is shown by the solid blue line, the linearcalculation by the red dashed line.with r MID = r = r =
9, the set of equations(36) was integrated forward in time. We chose the viscosityin these test simulations to be α = − ≪ h since we aremainly interested in the regime of propagating bending waves,and because the linear solution presented here becomes invalidfor large viscosities. The NIRVANA code was set up in an in-ertial frame using the expressions (10) and (18) to specify theinitial density and velocity, in combination with Eq. (17), with γ F ( r ) = δ ( r ) now being the disc tilt which is a function of r . Inthese simulations we used N r =
404 cell in radius, 5 cells perpressure scale height in the meridional direction (and di ff erentnumbers of vertical scale heights, depending on the initial warpamplitude, as described below) and N φ =
300 cells in azimuth.We used the same boundary condition as for the global simula-tions, which were described in section 3.1.
An example result is shown in Fig.1 which has a disc as-pect ratio of h = .
03 and maximum initial inclination angle δ MAX = ◦ . In this particular simulation the meridional domaincontains 15 pressure scale heights ( ∆ θ = . ◦ ) in total. Fig.1shows the inclination as a function of radius at di ff erent stagesof the evolution. The solid blue line represents the NIRVANAsolution, while the red dashed line shows the linear calculation.Moving from the top left to top right panel, we can observehow the initial pulse begins to broaden. At time t = .
96 wesee that the pulse has separated into an in-going and out-going bending wave. By time t = .
46 the pulses have reached theboundaries of the domain and have fully separated. It is clearthat NIRVANA is able to capture the evolution of these lowamplitude bending waves with a high level of accuracy.In Nelson & Papaloizou (1999) the numerical SPH and lin-ear solutions showed more substantial disagreement at the po-sition of the initial pulse ( r =
5) at the final stage of the simula-tion, with the separation of the two pulses not being captured asaccurately. This was probably due to the larger numerical di ff u-sion present in the SPH code, though it should be rememberedthat the SPH simulations were performed using only 20,000particles. The more accurate solution obtained with NIRVANA,on the other hand, is a clear indication that the numerical di ff u-sion is smaller.We also performed simulations for other values of the discthickness h and maximum initial inclination angle δ MAX . Thenumber of pressure scale heights used in the meridional direc-tion was scaled with the maximum inclination angle to avoidthe disc interacting with the meridional boundaries. The resultsare shown in Fig.2, and are presented at a stage of the simu-lations when the bending waves have propagated all the wayto the radial boundaries. This corresponds to the final panel ofFig.1. Since the bending waves propagate with half the soundspeed, the evolution occurs over a longer time in the thinnerdiscs ( h = . h = . .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 9 Fig.2 (top left and right panels), but the results are still in veryreasonable agreement for the higher inclination cases in Fig.2(lower left and right panels).We mention here that for even higher inclination anglesNIRVANA behaves more di ff usively, and the agreement be-tween the linear and non linear solutions gets worse. As foundby (Nelson & Papaloizou 1999), when the warps become non-linear they generate horizontal motions that are of the orderof the sound speed, creating shocks because of the symmetricform of the initial bending disturbance. The result is an e ff ec-tively higher viscosity operating in the disc, which is clearlynot accounted for in the linear calculations.
5. Tilted discs in binary star systems
We now discuss our results for the misaligned, tidally interact-ing discs. We set up the disc models according to the procedureoutlined in Sect. 3, and details of the model parameters can befound in Table 1. The models are characterized by the disc as-pect ratio, h , viscosity, α , and initial inclination with respect tothe binary orbit plane, γ F . The disc outer edge, R , was chosensuch that the disc is already very close to its tidal truncationradius, which is about one third of the distance between sec-ondary and primary star, D . In some models the initial radiusof the outer disc edge was varied also.In this study, we are interested in examining the evolu-tion of discs that fulfil the condition for wave-like warp prop-agation, h > α , as well as discs that support di ff usive warppropagation, h < α . For the models in the wave regime weset α = h , and for the runs in the di ff usive regime we set α = . ≫ h .The perturber is evolved on a circular orbit at a distanceof D =
30, and its mass is increased linearly from zero up to M S = M P over a time interval of 4 orbits, during which timeits distance is kept constant. Models 1 – 5 were performed inthe precessing frame. If we use Eq. (15) to estimate the preces-sion frequency, then we obtain a value of Ω F = − . · cos( γ F )[degrees / orbit]. We find, however, that a better fit to test simu-lations is given by Ω F = − . · cos ( γ F ) [degrees / orbit], whichis the value used in the simulations. The reason is that as thedisc gets tidally truncated it tends to precess slower, as can beseen from Eq. (15).The resolution was chosen such there are 5 cells per pres-sure scale height in the meridional direction for all models.Models 1 and 3 incorporated 15 scale heights in the meridionaldirection. Models 2, 4 and 5 used 22.5 scale heights (to allowfor a greater degree of twisting in the higher viscosity mod-els where warp communication is expected to be less e ffi cient).In the radial and azimuthal directions we used N r =
404 and N φ =
300 cells, respectively, for each of these models.Models 6 and 7 were performed in the non-precessing bi-nary frame. Since the disc midplane is inclined with respect tothe equatorial plane of the computational grid in these simu-lations, higher resolutions are necessary to avoid numerical er-rors. This is particularly true in the azimuthal direction becausetilting the disc causes a component of the disc vertical struc-ture to point along this direction. For a disc whose midplanecoincides with the equatorial plane of the computational grid, pairwise cancelation of fluxes ensures conservation of angularmomentum to machine accuracy. However, if the disc midplaneis inclined with respect to the equatorial plane of the computa-tional grid, the accuracy is limited by the advection scheme,and higher resolution in the azimuthal direction becomes nec-essary to avoid unphysical evolution of precession and inclina-tion angles. Hence we used 600 cells in azimuth. In order tobe able to resolve spiral waves for these extremely thin discs,we used 1056 cells in radius. The disc outer edge was chosento be a bit smaller than in the thick disc runs, since highly pro-nounced non linear e ff ects at the outer edge were seen in lowerresolution test simulations (described as Models 6a and 7a inTable 1). In e ff ect a strongly perturbed and narrow outer rim ofgas became partially detached from the main body of the discand developed a very di ff erent inclination and precession angleevolution in Model 6a, an e ff ect not observed in Model 6 witha smaller initial outer radius. This phenomenon is discussed inmore detail later in the paper. In order to accommodate an in-clined disc with an inclination of 10 ◦ we used 60 pressure scaleheights ( ∼ ± .
03 degrees) in the meridional direction, with 5cells per pressure scale height.We now describe the results of these models, where wehave ordered our discussion in terms of the disc aspect ratio, h . A column density plot for the disc in Model 1 is displayed inthe upper panels of Fig.3 corresponding to the end stage ofthis simulation. The left panel displays a projection of the disccolumn density onto the ( e , e ) plane of the binary referenceframe, and the right panel displays a projection onto the binary( e , e ) plane. At this stage the disc has precessed rigidly to β = −
180 degrees, so the disc appears edge on in the ( e , e )plane and face on in the ( e , e ) plane. The precession of 180degrees is evident, as the disc angular momentum vector nowhas a negative e component, whereas it had a positive e com-ponent in its initial state, as can be seen from the transformationgiven by Eq. (8) at t =
0. Another result of the interaction withthe secondary star is the formation of spiral waves, which areapparent in Fig.3.Fig.4 shows the evolution of the precession ( β ) and inclina-tion ( δ ) angles, respectively, for disc material orbiting at di ff er-ent radii. These have been calculated using Eqs. (22) and (21)respectively, where L has been calculated using Eq. (19), whichwe have integrated over radial shells of width ∆ r =
1. Thus thered dashed line in Fig.4 corresponds to disc gas between r = r =
2, the blue solid line to disc gas between r = r =
9, and the black dotted line to disc gas between r = r =
5. The solid black line shows the precession and inclina-tion angles of the entire disc. It is apparent that the lines are al-most indistinguishable in Fig.4 (upper left panel), showing thatthe di ff erent disc parts in Model 1 precess at the same uniformrate, such that the disc is in a state of rigid body precessionwith almost no twist. This behaviour is expected from lineartheory (Papaloizou & Terquem 1995). In Table 1 we show thewarp propagation timescale τ W and the di ff erential precession h α γ F R Frame τ W τ P Predicted Observedlabel behaviour behaviour1 0.05 0.025 45 9 precessing 5.43 48.01 rigid precession rigid precession2 0.05 0.1 45 9 precessing 10.87 47.60 rigid precession rigid precession3 0.03 0.015 45 9 precessing 9.05 43.90 rigid precession rigid precession4 0.03 0.1 45 9 precessing 30.19 47.02 rigid precession rigid precession5 0.03 0.1 25 9 precessing 30.19 36.37 rigid precession rigid precession6 0.01 0.005 10 8 binary 22.77 33.87 rigid precession rigid precession7 0.01 0.1 10 8 binary 227.7 36.25 disrupted / broken rigid precession6a 0.01 0.005 10 10 binary 31.83 21.48 disrupted / broken broken7a 0.01 0.1 10 10 binary 318.32 19.12 disrupted / broken disrupted Table 1.
Table of runs. The disc is characterized by the aspect ratio, h , viscosity parameter, α , and initial inclination angle γ F .We also list the initial disc outer radius R and the reference frame in which the simulation was performed. The warp propagationtimes, τ W and di ff erential precession times, τ P , for each of the models are also listed and indicate the theoretically expected andobserved behaviours. Fig. 3.
Column density plot at time t = . e , e ) binary plane, and the right panels to projection onto the ( e , e ) plane. .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 11 Fig. 4.
Evolution of precession (left panels) and inclination angles (right panels) for Models 1 and 2. The blue solid line corre-sponds to disc material between r = r =
9, the red dashed line to disc material between r = r =
2, and the blackdotted line to disc material between r = r =
5. The black solid line represents the entire disc.timescale, τ P . In the bending wave propagation regime ( h > α ), τ W = R / c S , since warps propagate with half the sound speedin this regime. In the di ff usive regime ( h < α ) τ W = α R /ν .The warp propagation timescale should be compared to the dif-ferential precession timescale τ P ∼ / Ω F . If τ W < τ P it isexpected that the disc will precess as a rigid body, while if τ W > τ P the disc may become disrupted as a result of di ff er-ential precession. Hence, we see that the rigid precession ofthe disc in Model 1 agrees with expectations. Our Model 1 isvery similar to Model 1 in Larwood et al. (1996), which alsoshows rigid body precession. The inferred precession period inour Model 1 is about 146 orbits, which is a little larger than theperiod given by Eq. (15). The di ff erence arises because the discis tidally truncated and shrinks slightly, so its precession ratedecreases.Looking at the upper right panel of Fig.4, we observe thatthe inclination for Model 1 changes only slightly during thesimulation. The oscillation seen in the inclination rates is drivenby the secondary star. Since the orbital plane is inclined with re-spect to the disc midplane, the vertical component of the grav-itational force due to the secondary star causes the disc to per-form a forced oscillation with a frequency equal to twice the bi-nary frequency. Papaloizou & Terquem (1995) argued that thedisc is expected to align with the binary orbit on a time scalesimilar to that over which binary torques cause the total angu-lar momentum of the disc to evolve. If the viscously-inducedoutward expansion of the disc is balanced by tides due to the companion, then disc-alignmnent should occur on the viscousevolution time, T ν . For Model 1 we estimate T ν = R /ν ∼ ∼ ǫ , which is a measure of the disc thickness, and isroughly equivalent to √ H / R . Their model a displayed in theirFigure 7 has H / R ≃ . α = .
01 and D / R = .
3, similar toour model 1. Although a detailed comparison is not possible,it appears that there is general agreement between linear the-ory and our non-linear simulation, since both show almost nodistortion of the disc due to twisting or warping.Column density plots for Model 2 are displayed in thelower panels of Fig.3. Warp propagation in this higher viscos-ity model ( α = .
1) is expected to be di ff usive, and thereforeless e ffi cient than for Model 1 (as shown by the warp propa-gation times listed in Table 1). This run was also evolved untilthe precession angle reached -180 degrees. The high viscosityoperating in the disc leads to the damping of the spiral wavesbefore they can propagate very far, and thus they are not appar-ent in the images.The precession angles for Model 2 are displayed in Fig.4(lower left panel). We observe that this more viscous disc de-velops a small di ff erential twist before it attains a state of Fig. 5.
Column density plots for Model 3 (upper panels) and Model 4 (middle panels) at time t = .
0, and Model 5 at time t = .
0. The left panels correspond to projection onto the ( e , e ) plane, the right panel to projection onto the ( e , e ) plane.rigid precession. This suggests that the disc needs to developa slightly more distorted shape in order for internal stresses tocounterbalance the companion-induced di ff erential precessionwhen warp communication is less e ffi cient. Nonetheless, rigidbody precession is expected in this case, and the numerical re-sults are in agreement with this expectation. Examining the inclination evolution for Model 2 shown inthe lower right panel of Fig. 4, we see it is very similar tothat observed for Model 1. The viscous time scale in Model2 should be approximately 4 times shorter than for Model 1,but there is no indication that the inclination evolution is fasterin this case. We note that the inclination evolution rates do notappear to have reached final steady values for either Models 1 .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 13 or 2 by the end of the simulations, such that an accurate assess-ment of the longer term evolution rate is not possible. This maybe because the outer disc radius has not had time to equilibrateduring the simulations as they have only run for a fraction of theglobal viscous timescale of the disc. We are able to conclude,however, that the inclination evolution occurs over timescalesmuch longer than the precession time, and appears to be con-sistent with evolution on the global viscous time of the disc. As shown in Table 1, Models 3, 4 and 5 all have h = . α = . via bending waves, whereas Models 4 and 5 have α = . ff usively. Table 1 also shows that the warppropagation for Models 4 and 5 is expected to be considerablyless e ffi cient than for Model 3. Models 3 and 4 had inclinationsof 45 ◦ , whereas Model 5 had an inclination of 25 ◦ , which in-duces a more rapid precession of the disc, and therefore has agreater tendency to twist the disc up.At the outset of this project it was our intention to run allmodels until they had precessed by 180 ◦ . For models 3 and4, however, we found that the discs become eccentric on atimescale of approximately 60 orbits (16 binary orbits), lead-ing to undesirable changes in the disc structure due to our useof an open inner boundary condition. Thus we did not evolvethese models until they had precessed by 180 ◦ . We note that thetimescale for the observed disc eccentricity growth correspondsclosely to that obtained by (Kley, Papaloizou & Ogilvie 2008).Model 5 did not su ff er from this problem because the lowerinclination angle induced more rapid precession, such that thismodel could be evolved until it had precessed by 180 ◦ , but priorto growth of significant eccentricity. The study of disc eccen-tricity goes beyond the scope of this paper, but the results formodels 3 and 4 suggest that very long term integrations maylead to the formation of inclined and eccentric discs. The factthat we did not evolve models 3 and 4 through 180 ◦ of preces-sion means that the upper and middle panels of Fig.5 displayimages which correspond to 90 ◦ of precession.It is clear from these panels that the disc in Model 3has precessed as a rigid body, with almost no twist apparent.Consequently, after precessing through 90 ◦ the disc appearsedge on when projected into the ( e , e ) plane and face onwhen projected into the ( e , e ) plane. The middle right panelof Fig. 5, however, shows that the disc in Model 4 has devel-oped a small twist due to the inner disc not having precessedby 90 ◦ at the same moment when the outer disc has. The lowerpanels of Fig. 5 show column density plots for Model 5, whichwas evolved until the outer disc had precessed by 180 ◦ (assistedby the fact that the precession is faster in this case). It is appar-ent from these panels that the disc has developed a small twistsince the inner disc has not precessed a full 180 ◦ by the end ofthe simulation.The precession angles for Model 3 are plotted in the top leftpanel of Fig. 6. The blue solid line corresponds to disc materialorbiting between radii r = r = r = ff ectively lie on top of each other shows that the disc precessesas a rigid body with essentially no twist. Rigid-body precessionis expected in this case since the warp propagation time is lessthan the di ff erential precession time (see Table 1).The inclination angles for Model 3 are plotted in the topright panel of Fig. 6. The line styles and colours correspondto the same disc annuli as described above. It is clear that theinclination evolution in this case is very slow. The global vis-cous evolution time for this disc is approximately 10,000 orbits,and a linear extrapolation of the inclination evolution shown inFig. 6 gives a time of ∼ ff erential precession, with the outer discprecessing faster than the inner disc. The magnitude of the twistis approximately 12 ◦ . Thereafter the disc is seen to precess asa rigid body, maintaining the same twisted shape. We note thatthis disc is expected to achieve rigid body precession from con-sideration of the warp propagation and di ff erential precessiontimes (see Table 1).The inclination angles for Model 4 are shown in the rightmiddle panel of Fig. 6. We see that the disc in this case has de-veloped a small warp, with the inner and outer discs having adi ff erence of inclination of about 1 ◦ . We also observe that theinclination is evolving more rapidly than is the case for Model3. The global viscous evolution time for this disc is approxi-mately 1500 orbits, and extrapolating the inclination evolutionobserved in Fig. 6 gives an estimated time for alignment of ∼ h = . D / R = ffi cult to know the pre-cise value of the α viscosity operating in the SPH simulations,it is likely that the e ff ective α value lies between the valuesused in our Model 3 ( α = . α = . ff erent behaviour, with ourmodels quickly achieving a state of rigid body precession, witha small twist and warp which vary smoothly across the disc.Model 9 of Larwood et al. (1996) shows significant disruptionof the disc, which e ff ective breaks discontinuously into two dis-tincts parts which become mutually inclined by approximately25 ◦ . The origin of this discrepancy is unclear, but one possi-bility is that SPH simulations using 20,000 particles are onlymarginally able to resolve the vertical structure of a disc with h = .
03. This may have the e ff ect of reducing the e ffi cacy ofwarp propagation below that which is expected from consider-ation of the disc physical parameters.We now turn to Model 5. The precession angles for thisrun are displayed in the lower left panel of Fig.6. Given thatthe precession frequency, given by Eq. (15), is proportional Fig. 6.
Evolution of precession (left panels) and inclination angles (right panels). The blue solid line corresponds to disc materialbetween r = r =
9, the red dashed line to disc material between r = r = r = r =
5. The solid black line represents the entire disc. In the lower left panel the free particle precessionrates for the inner (dashed-dotted line) and outer disc (dashed-triple-dotted line) are also depicted for Model 5. These have beencalculated using Eq. (38).to cos ( γ F ), the disc precesses faster in this simulation than inModels 1 – 4, and this faster precession can be observed inFig. 6. The lower left panel of this figure shows the preces-sion angles as a function of time for disc annuli located be-tween r = r = r = ff erential precession inthe disc, and the consequence of this is that Model 5 showsa slightly larger twist (15 ◦ ) than Model 4 (12 ◦ ). Also plottedin this panel are precession angles that would be observed ifthe inner (black dashed-dotted line) and outer (black dashed-triple-dotted line) disc annuli precessed at their free particle rates given by (Papaloizou & Terquem 1995): ω Z = − Ω K GM S D (3 cos δ − . (38)As can be seen from the figure, the disc parts tend initiallytowards their free particle rates, such that the outer disc pre-cesses faster than the inner disc. Information about this di ff er-ential precession is communicated across the disc, and internalstresses are established which cause the inner disc precession tospeed up and the outer disc precession to slow down, such thatprecession at a uniform rate occurs after approximately 10 or-bits. The warp communication and di ff erential precession timesgiven in Table 1 for this model predict rigid-body precession, as .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 15 Fig. 7.
Precession rates due to the contributions discussed in Sect. 3.6 for Model 5. The blue solid line corresponds to disc materialbetween r = r =
9, the red dashed line to disc material between r = r = r = r =
5. These rates have been calculated using Eqs. (34). In the lower right panel the total rate is displayed,which has beeen calculated using Eq. (33), where we subtracted the constant mean precession rate Ω F of the precessing frame. Fig. 8.
Inclination rates due to the contributions discussed in Sect. 3.6 for Model 5. The blue solid line corresponds to discmaterial between r = r =
9, the red dashed line to disc material between r = r = r = r =
5. These rates have been calculated using Eqs. (34).observed. The emergence of internal stresses as the disc twistsup is illustrated in Fig.7, where the di ff erent physical contribu-tions to the precession rate are displayed for Model 5. Thesehave been calculated according to the procedure outlined inSect. 3.6. Here the di ff erent linestyles correspond to the samedisc annuli described above (solid blue line: r = r = r = Ω F < Fig. 9.
Time averaged inclination change rates as function ofradius for Model 5. Red dashed-dotted line: radial flux of angu-lar momentum; blue dotted line: gravitational interaction withthe companion star; green dashed line: viscous friction betweenadjacent radial disc shells; black dashed-triple-dotted line: totalinclination change rate.outer annulus would lead the precession of the reference frame,if the companion’s gravity was the only contributor to the to-tal precession rate. In fact, this is what is observed until about t =
10 orbits, during which time the disc precesses di ff eren-tially (see Fig.6, lower left panel). As the disc starts to developa twist due to di ff erential precession, misalignment of disc mid-planes as a function of radius causes pressure forces to generatehorizontal shear motions that drive the propagation of bendingdisturbances, and these communicate information about pre-cession angle changes across the disc. The resulting internalstresses are clearly proportional to the degree of twisting, al-lowing the disc to achieve a state of rigid body precession oncea su ffi ciently twisted disc has been established. At this point thedisc structure becomes quasi-static in the precessing frame aseach disc annulus precesses at the same rate. This sequence ofevents is demonstrated by the upper left panel of Fig.7, whichshows the contribution to the precession rate from the pressure-induced radial flux of angular momentum. After t =
10 orbits,the contribution to the precession rate from the radial flux al-most completely compensates the e ff ect of gravity, such thatthe total precession rate (relative to the precession of the ref-erence frame) is zero for all disc annuli, as can be seen in thelower right panel of Fig.7, which shows the sum of the contri-butions to the precession rate. This demonstrates that the discmanages to redistribute the angular momentum such that rigidbody precession is achieved. Once again we note that in thisfigure we have subtracted the mean precession rate Ω F of theprecessing frame.The e ff ect of viscosity as calculated in Sect. 3.6 is depictedin the upper right panel of Fig.7. We observe that the e ff ectof viscous friction between di ff erentially precessing radiallyadjacent disc shells causes the disc to precess more rigidly.However, the dominant e ff ect of the viscosity is the damping ofthe horizontal shear motions induced by the twist. Since these shear motions are responsible for driving the bending distur-bances, the viscosity will lead to weakened communication,and the redistribution of angular momentum due to radial ad-vective fluxes shown in the upper left panel of Fig.7 will be sup-pressed. This e ff ect is dominant and opposite to the one seen inthe upper right panel of Fig.7.The inclination angles for Model 5 are shown in the lowerright panel of Fig. 6. It is clear that the disc develops a smallwarp, with the inner disc having ∼ ◦ higher inclination thanthe outer disc. The rate of global inclination change is similarto that observed in Model 4, suggesting alignment of the discin the global viscous evolution timescale.We plot the di ff erent physical contributions to the rate ofinclination change for Model 5 in Fig.8. In Fig.9 we also dis-play the time averaged contributions to the rate of inclinationchange as a function of radius, where the time averaging wasperformed over the full length of the simulation. As may be ob-served in the lower left panel of Fig.8, and even more clearlyin Fig.9 (blue dotted line), the companion’s gravity induces aglobally negative net inclination change, leading to coplanarityof the entire disc on the viscous time scale, in agreement withthe lower right panel of Fig.6. This is consistent with the find-ings of Papaloizou & Terquem (1995) that the zero-frequency(secular) gravitational interaction will lead to coplanarity. Wenote that because of the time averaging, this is the dominantgravitational contribution that we pick up in Fig.9 (blue dottedline).More interestingly, we can observe from the upper leftpanel of Fig.8, and in Fig.9 (red dashed-dotted line), that thecontribution from the radial fluxes not only counteracts the ef-fect of gravity, it actually overshoots slightly, causing the innerdisc to become more inclined than the outer disc. This changeof inclination is seen in the upper right panel of Fig.6. Thus thecommunication of bending disturbances on the one hand leadsto rigid precession of the disc, and on the other hand causes thedisc to become slightly warped. Viscosity tends to reduce theamount of warping, such that it acts to flatten the di ff erentialinclination profile across the disc radius, as can be seen fromthe green dashed line in Fig.9. Thus the total average inclina-tion rate in Fig.9 (dashed-triple-dotted line) is negative for theentire disc, while the magnitude is slightly bigger for the outerdisc, with the consequence that the disc has developed a slightwarp and tends to become coplanar with the orbital plane onthe viscous timescale, as seen in the lower right panel of Fig.6. We now discuss Models 6 and 7, for which the discs have as-pect ratio h = .
01, an inclination angle of 10 ◦ , and were per-formed in a non precessing frame. Model 6 has α = . via bending waves in this case, and Model7 had α = .
1, so warps propagate di ff usively in this model.We see from Table 1 that di ff erential precession is expectedto disrupt the discs in Model 7, as the warp propagation timeis longer than the di ff erential precession time. Rigid-body pre-cession is expected for Model 6. .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 17 Fig. 10.
Column density plots for Model 6. The left-hand panels are projections onto the ( e , e ) plane, and the right-hand panelsare projections onto the ( e , e ) plane. The disc structure is displayed at the di ff erent times indicated. Fig. 11.
Evolution of the precession (left panels) and inclination angles (right panels). The blue solid line corresponds to discmaterial between r = r =
8, the red dashed line to disc material between r = r =
2, and the black dotted line to discmaterial between r = r =
5. The solid black line represents the entire disc.Column density plots showing the disc in Model 6 are pre-sented in Fig.10 for di ff erent times during the evolution. Thetop panels show the initial state of the disc. The disc appearsedge on in the ( e , e ) plane and face on in the ( e , e ) plane.As we move down to the next panels the outer edge of the dischas precessed by 45 degrees at time t = .
4, as can be seenin Fig.11 (upper left panel). The inner region of the disc, how-ever, has only precessed by about 12 degrees at this stage andso appears almost edge-on in the left panel. In the third panelsdown the outer disc has precessed by 90 degrees at t = . e , e ) plane and face on in the ( e , e ) plane. The inner disc has precessed by approximately 70degrees by this time, and so does not appear edge-on in eitherprojections. The bottom panels show the disc structure at thefinal output of the simulation which was halted at t =
44 orbits,by which time the outer disc had precessed by approximately160 ◦ .The precession angles for Model 6 are shown in the upperleft panel of Fig. 11. The blue solid line corresponds to discmaterial orbiting in an annulus between r = r = r =
1– 2. As with Models 4 and 5, we see that the disc inner andouter annuli begin to precess at the local free-particle rates, andthe disc develops a significant twist which reaches a maximumamplitude of approximately 30 ◦ . As the disc becomes increas-ingly twisted, however, internal stresses are established which cause the inner disc precession rate to increase and the outerdisc precession rate to decrease. Eventually the disc reaches astate of rigid body precession after a time of about 10 orbits.After this time we see that there is a long term readjustment ofthe degree of twist such that by t =
44 orbits the twist angle hasreduced from ∼ ◦ to ∼ ◦ .In Fig. 12 we display the time integrated contributions tothe precession rate as a function of radius for Model 6, follow-ing the procedure described in Sect. 3.6. Note that the average(negative) precession rate of the disc has been subtracted inthis figure, such that a positive value corresponds to preces-sion that will lag that of the overall disc, and a negative valuecorresponds to precession that leads. The contribution from thecompanion’s gravity is shown by the blue dotted line, and thisis very similar to the free-particle precession rate which is plot-ted using the black solid line. The red dashed-dotted line showsthe contribution from the pressure-induced radial flux of an-gular momentum, and this is seen to almost exactly counter-balance the e ff ect of gravity, showing that it is the pressure-induced radial motion in the disc which is responsible for halt-ing the di ff erential precession and inducing rigid body preces-sion. The contribution from the disc viscosity is shown by thegreen dashed line, and is seen to be negligible.The inclination angles for Model 6 are shown in the topright panel of Fig. 11, where the line styles and colours cor-respond to the disc annuli described above. As with Model 3,we see that the disc develops a warp such that the inner disc .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 19 Fig. 12.
Average precession rates as function of radius forModel 6. Red dashed-dotted line: radial flux of angular mo-mentum; blue dotted line: gravitational interaction with thecompanion star; green dashed line: viscous friction betweenadjacent radial disc shells; black dashed-triple-dotted line: to-tal inclination change rate. The black solid line represents theprecession rate of free particles.
Fig. 13.
Average inclination change rates as function of radiusfor Model 6. Red-dashed-dotted line: radial flux of angularmomentum; blue dotted line: gravitational interaction with thecompanion star; green dashed line: viscous friction between ad-jacent radial disc shells black dashed-triple-dotted line: totalinclination change rate.has a larger inclination to the binary plane than the outer disc.Indeed, the inclination of the inner disc is seen to increase ini-tially, before decreasing after approximately 20 orbits as thedisc as a whole aligns with the binary orbit plane. As discussedpreviously, disc alignment is expected to occur on the viscoustime scale, which for Model 6 is approximately 3 × orbits.It is obvious from Fig. 11, however, that the disc in our simula-tion is not aligning on such a long time scale, but will actuallyalign within ∼
650 orbits. Possible reasons for this enhancedalignment rate are discussed below. Fig. 13 shows the time integrated contributions to the rateof inclination evolution as a function of radius for Model 6.We see that the companion’s gravity causes the inclination todecrease in the disc inner and outer regions, and in the inner re-gions the radial flux of angular momentum opposes this and in-creases the inclination during the simulation. Indeed the radialflux contribution causes the inclination angle to grow there, asobserved in the lower right panel of Fig. 11. The outer parts ofthe disc experience negative contributions from both the com-panion’s gravity and the pressure-induced radial flux, drivingalignment of the disc on the relatively short timescale describedabove. It is clear that the disc viscosity plays a negligible rolein driving the inclination evolution (green dashed line).We may compare our model 6 with model c in Figure 7 ofLubow & Ogilvie (2000), which had H / R ≃ . α = . D / R = .
3. Detailed comparison between the results pre-sented in Lubow & Ogilvie (2000) and our model 6 is not pos-sible, but we note that there is general agreement in the pre-diction that a disc with H / R ≃ .
01 and α . H / R will showsignificant distortion due to twisting and warping. A detailedcomparison between our non linear simulation results and thepredictions of linear theory will be presented in a future publi-cation.We now discuss the results from Model 7. In Fig.14 we dis-play column density plots for this model. The top panels showthe initial state. The second panels down show the disc at time t = .
5. At this stage the outer disc edge has precessed by 90degrees and appears edge on in the ( e , e ) plane, and face onin the ( e , e ) plane. By contrast, the inner disc has not pre-cessed very much and still appears almost edge on in the ( e , e ) plane and face on in the ( e , e ) plane. At time t =
44 orbitsthe outer disc has precessed by 180 degrees and appears edgeon in the ( e , e ) plane and face on in the ( e , e ) plane. Theinner disc, however, has only precessed by about 70 degreesby this time, and so does not appear edge on in either projec-tion. The bottom panels show the final stage of the simulationwhen the outer disc has precessed by about 220 degrees. It isclear from these images that the disc in Model 7 has adopted ahighly twisted shape, but one in which the precession angle fordi ff erent disc annuli varies very smoothly across the disc.The precession angles for Model 7 are displayed in thelower left panel of Fig.11. We observe that the disc is in a stateof di ff erential precession at the end of the simulation, and thetwist between inner and outer disc annuli is approximately 110degrees at this stage. Although this twist is large, the actual rateof di ff erential precession at the end of the simulation is verysmall, as the precession rates of the inner and outer disc haveconverged during the run. Indeed, extrapolation of the preces-sion angles in Fig. 11 indicates that rigid body precession willbe achieved at a time equal to t = . . ◦ . This suggests that even thoughthe warp propagation time is much smaller than the di ff erentialprecession time for Model 7 (see Table 1), internal stresses canbe set up within the disc that cause uniform precession oncethe disc becomes highly twisted. Adopting a definition of discdisruption which requires the twist to become greater than 180degrees, we suggest that the disc in Model 7 will not be dis- Fig. 14.
Column density plots for Model 7. The left-hand panels are projections onto the ( e , e ) plane, and the right-hand panelsare projections onto the ( e , e ) plane. The disc structure is displayed at the di ff erent times indicated. .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 21 Fig. 15.
Diagram illustrating the alignment torques which arisein highly twisted discs, as discussed in the text.rupted, but will instead adopt a highly twisted but uniformlyprecessing configuration.The inclination angles for Model 7 are shown in the lowerright panel of Fig. 11, where the lines styles correspond to thedisc annuli described above. We can see that the disc in thiscase develops a small warp with the di ff erence between the in-clinations of the inner and outer disc being ∼ ◦ . We also notethat it appears that the integrated tilt of the whole disc is ac-tually smaller than for any of the individual disc annuli, butthis is an e ff ect of averging over a significantly twisted disc. Asobserved for Model 6, we see that the inclination of the discevolves rather quickly. The global viscous evolution time forModel 7 is ∼ ∼
100 orbits.What is the explanation for this rapid alignment observedfor Models 6 and 7 ? One possibility that we have explored isthat numerical di ff usion may cause a tilted disc to align withthe equatorial plane of the computational grid, since the advec-tion routine does not conserve total angular momentum whenthe disc azimuthal velocity is not directed along the azimuthaldirection of the computational grid. Lower resolution test cal-culations performed for a tilted disc without a companion indi-cate that this e ff ect can cause eventual disc alignment, but at arate which is at most 7 times smaller than observed in Models6 and 7. This suggests that numerical di ff usion is not the causeof the rapid alignment.Another possibilty arises when we consider the evolutionof a rigidly precessing, tilted disc which is in a precessing ref-erence frame whose precession frequency is equal to that of thedisc. Consider first the situation in which the angular momen-tum vectors for all disc annuli lie in the y − z plane of the binaryframe initially (i.e. the z direction of the precessing frame liesinitially in the y − z plane of the binary frame, pointing alongthe positive y and z directions). A torque exerted on disc annuliin the y direction, (a ‘ y -torque), will cause the disc inclinationto change. A torque exerted in the x direction (an ‘ x -torque’),will cause precession at a rate which di ff ers from the preces-sion rate of the frame. A positive x -torque causes precessionto occur at a rate which is faster than the frame, and a nega-tive x -torque causes slower precession. This general picture is expressed mathematically by Eqs. (34), and is shown diagra-matically in the left hand image of Fig. 15. For a disc whichprecesses uniformly without any twist, then the companion’sgravity will exert positive x -torques on the outer disc annuli,and negative x -torques on the inner disc annuli, as it tries tocause the disc to precess di ff erentially. Internal stresses, how-ever, will exactly counterbalance these x -torques to maintainuniform precession. In this simple scenario there is no torqueexerted on the disc in the y direction to modify its inclination.Consider now the case of a highly twisted disc which is be-ing forced to precess uniformly, and whose annuli all have thesame inclination. Gravity will again try to induce di ff erentialprecession, but now the torque exerted on disc annuli will haveboth an x and a y component due to the twist. Resolving thevectors associated with the internal stresses required to enforceuniform precession shows that their combined x -torques oper-ate in opposite directions such that they cancel when integratedover the disc. The combined y -torques, however, operate in thesame direction, implying that the internal torques generate a net y -torque on the disc. Conservation of angular momentum ob-viously prohibits this from arising, suggesting that a y -torquemust operate, whose e ff ect is to change the disc inclination.This is shown diagramatically in the right panel of Fig. 15.Given that the internal torques are operating on the precessiontimescale in order to maintain uniform precession, this argu-ment suggests that a highly twisted disc may align on the pre-cession time, as is observed for Models 6 and 7. Interestingly,the resultant y -torque that arises in this picture is proportionalto sin ( ∆ β ), where ∆ β is the twist amplitude, and the rate of in-clination evolution observed in Models 6 and 7 is in agreementwith this scaling.If the above argument is correct, then it suggests that discswith modest levels of twist will align on the global viscoustimescale of the disc, but highly twisted discs may align onthe much shorter precession timescale. This obviously has sig-nificant implications for astrophysical systems which may beexpected to harbour highly twisted discs, such as the X-ray bi-naries Her X-1 and SS433. We now briefly discuss the results of two lower resolution mod-els which were run during an early stage of this project: Models6a and 7a. These were identical to Models 6 and 7 except thatthe number of grid cells in the radial direction was N r = R =
10 instead of R = ff ect of in-creasing the rate of di ff erential precession across the disc, pro-vided that tidal truncation does not cause the disc to shrink backdown to R = t = . r ∼ r ∼
10, has e ff ectively de-tached itself from the main body of the disc, and has evolved Fig. 16.
Column density plots of the two lower resolution runs Models 6a and 7a. The upper panels show two di ff erent projectionsfor Model 6a at time t = . ff erent projections for Model 7a at time t = . ff erent precession and inclination anglescompared to the main body of the disc.The inclination and precession angles for Model 6a are dis-played in the upper left and right panels of Fig. 17, respectively.The precession angle of the annulus which lies between radii r = r =
7. The solid blue line corresponds tothe disc annulus lying between r = r = r = r = r =
8, with thedetached outer rim having precessed significantly faster thanthe inner disc. The inclination angles plotted in the upper rightpanel show that the disc becomes significantly warped, with the outer rim tending to align with the binary and the innerdisc actually increasing its inclination. Although the change ininclination appears to occur fairly smoothly across the disc, itshould be noted that the disc has developed a fairly large warpof 5 ◦ .The behaviour of Model 6a is reminiscent of the brokendisc presented by Larwood et al. (1996) (their broken Model 9had h = .
03, inclination of 45 ◦ , D / R =
3) except that the dis-continuous break in Model 6a occurs near to the edge of thedisc. It would appear that Model 6a breaks because the largerdisc radius induces a larger rate of di ff erential precession, andalso possibly because warp propagation is retarded in the outerdisc regions which are subject to non linear density structuresand shocks because of strong interaction with the binary com-panion. .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 23 Fig. 17.
Evolution of precession (left panels) and inclination angles (right panels) for Models 6a and 7a. The blue solid linescorrespond to disc annuli of unit radius lying between r = r =
10. The solid blue line denotes the annulus lying between r = r =
10. The red lines correspond to disc annuli of unit width lying between r = =
7. The solid red line denotesthe annulus lying between r = r = t = ∼ ◦ ), butdoes not show any of the discontinuous behaviour observedin Model 6a. Instead the twist varies very smoothly across thedisc. The precession and inclination angles for Model 7a areshown in the lower left and right panels of Fig. 17, respectively.The solid red line shows the precession angle of the annulus lo-cated between r = r = ∼ ◦ , and the precession rates of innerand outer disc remain very di ff erent such that this disc will dif-ferentially precess through 180 ◦ . So, by our working definitionof disc disruption, this disc will be disrupted because of di ff er-ential precession. The inclination angles show that the disc alsobecomes significantly warped, with the outer disc aligning withthe binary orbit more rapidly than the inner disc.To summarise: Model 6a and 7a both show examples ofdiscs which break or are disrupted because of di ff erential pre-cession. In the case where warps propagate via bending waves,the disc breaks discontinuously into two independently pre-cessing parts, but which themselves remain as coherent struc-tures which precess as rigid bodies. In the case where warpspropagate di ff usively, viscosity is able to smooth out any dis-continuous behaviour, and instead the disc becomes disruptedbecause the twist exceeds 180 degrees. In all probability thedisc annuli in Model 7a will continue to precess di ff erentially inthe longer term, substantially destroying any coherent disc-likestructure. The larger disc radius in this case, when comparedwith Model 7, causes the di ff erential precession rate across thedisc to be larger, and thus explains why Model 7a shows dis- ruption, whereas Model 7 results in a highly twisted disc whichis able to attain a state of rigid-body precession.
6. Conclusions
We have presented non linear simulations of accretion discsin close binary systems in which the binary midplane is mis-aligned from the binary orbital plane. Previous work on thisproblem, which employed linear theory, semi-analytic tech-niques and SPH simulations (Papaloizou & Terquem 1995;Larwood et al. 1996; Lubow & Ogilvie 2000), suggests thatunder suitable conditions the disc should become mildlywarped and precess as a rigid-body around the angular momen-tum vector of the binary system. The main aim of the presentstudy is to examine in detail the structure and evolution ofdiscs in misaligned binary systems, subject to variation of theimportant physical parameters. We used the grid-based codeNIRVANA to perform the numerical simulations, and in prin-ciple this should have the advantage over previous SPH simu-lation studies that the disc viscosity is a well defined and con-trollable parameter.Following Nelson & Papaloizou (1999), the propagation oflinear bending waves was used as a means to calibrate the codeand test its ability to propagate warps. Direct comparison withlinear calculations show that the code can propagate bendingwaves with a high degree of accuracy for a range of disc modelsand warp perturbation amplitudes.In our study of discs in misaligned binary systems, a num-ber of di ff erent models were considered. Discs with aspect ra-tios h = .
05, 0.03 and 0.01 were studied, and for each ofthese values models were computed with viscosity parameters α = h / α = .
1. The smaller value of α allows bending waves to propagate, whereas the larger value causes warps toevolve di ff usively.It is expected that discs for which the warp propagationtime, τ W , is shorter than the di ff erential precession time, τ P ,will undergo rigid-body precession. Discs which do not satisfythis criterion may be disrupted by strong di ff erential preces-sion. All our models with h = .
05 and 0.03 have τ W < τ P , andthe simulations show that each of these discs attains a state ofrigid-body precession. Discs with lower viscosity maintained aprecessing structure with almost no discernible twist and warp,due to the highly e ffi cient warp propagation associated withbending waves. Discs with larger viscosity propagate warpsless e ffi ciently, and these models underwent a short period ofdi ff erential precession, setting up a twist in the disc, prior toattaining a state of rigid-body precession. Our analysis showsthat internal stresses are established in the disc as it becomestwisted, which transport angular momentum across its radiusand cause it to precess rigidly. Models in which warp propaga-tion is the least e ffi cient develop the largest twist, since the in-ternal stresses which hold the disc together against di ff erentialprecession are proportional to the twist amplitude. Examinationof the simulations shows that the dominant contributor to theseinternal stresses are pressure-induced radial angular momen-tum fluxes which are driven by local misalignment of disc mid-planes. When α = .
1, the final amplitude of the twist was justa few degrees for the h = .
05 disc, and for the h = .
03 disc itwas between 12 – 15 degrees, depending on the binary inclina-tion. The low viscosity models showed essentially no twist.In addition to becoming twisted, a disc may also becomewarped. We find that for all models with h = .
05 and h = . ff erences in in-clination across all disc annuli being less than 1 degree. It isexpected that these discs will also align with the binary or-bital plane on the global viscous evolution time of the disc(Papaloizou & Terquem 1995; Larwood et al. 1996), and allsimulations with h = .
05 and 0 .
03 are fully consistent withthis expectation.We also considered a number of thin disc models with h = .
01, where we not only varied the viscosity but also theouter radius (taking values of either R = R = R = ◦ before achieving a state of solid body pre-cession, as expected since τ W < τ P for this model. A modelrun with α = . R = ff erential precession since τ W > τ P in this case. Duringthe simulation it developed a very large twist ∼ ◦ , but as thedisc become increasingly twisted the rate of di ff erential preces-sion decreased. At the end of the run the disc was experiencingvery slow di ff erential precession, and extrapolation of the rateof twisting indicates that this disc will achieve rigid-body pre-cession with a twist of ∼ ◦ . This is the first indication inour results that a disc which is predicted to be disrupted cannonetheless generate internal stresses to enforce rigid preces-sion when the degree of twist becomes large.Models with h = .
01 and R =
10 showed di ff erent be-haviour, largely as a result of the larger precession rate whichis induced for a larger disc, and suggest that the models runwith R = α = h /
2, which supportsbending waves, was observed to undergo modest tidal trunca-tion and to break into two distinct disc parts which precessedindependently of one another, an outcome which is similar toone obtained using an SPH simulation by Larwood et al. (1996)but for a disc with somewhat di ff erent parameters. A narrowrim at the edge of the disc, running between radii r = R = τ W > τ P for this larger disc model, such that the breaking ofthe disc is indeed expected.The model with α = . ff usively, andwas observed to di ff erentially precess through 180 degrees, butmaintained a very smooth twist profile throughout the simula-tion. This disc is clearly disrupted through strong di ff erentialprecession, and the internal stresses are insu ffi cient for the discto achieve rigid-body precession. Over longer times this discshould become completely twisted up by di ff erential preces-sion such that it loses all semblance of a disc-like structure.The results for the R =
10 discs are interesting, as theyindicate two di ff erent modes by which a disc may break or bedisrupted. In the bending wave regime, our results suggest thatdisc disruption may occur via a discontinuous change in theprecession rate at a particular radius, such that the disc breaksinto two distinct pieces which precess independently. Withineach separate disc-piece warp communication allows for rigid-body precession. The detachment of the disc outer rim may inpart be influenced by the strong tidal interaction there, sincenon linear features and shocks may have the e ff ect of reducingthe e ffi cacy of warp communication. In the di ff usion regime, itappears that viscosity is able to smooth out any discontinuities,and the disc is disrupted through global di ff erential precessionof the disc annuli, with a smooth twist profile being maintainedacross the disc radius.A significant result obtained for the h = .
01 modelsis that these highly twisted discs align with the binary orbitplane much more quickly than the thicker discs, which alignon the viscous timescale. Although we find that numerical dif-fusion can cause an inclined disc to align with the equatorialplane of the computational domain, this occurs on much longertimescales than are observed. We tentatively suggest that align-ment of highly twisted discs occurs on the precession timerather than the viscous time, causing rapid alignment of verythin discs.We have not considered di ff erent binary mass ratios, insteadwe restricted our study to the case of an equal mass binary.As Eqs. (38) and (15) indicate, the free and mean precessionrates scale linearly with the companion mass for a disc of fixedouter radius. Hence we would expect the di ff erential precessiontimescale τ P to be increased for less massive binary compan-ions, and consequently the di ff erential twist observed in ourmodels should be reduced in such systems.Our work has a number of astrophysical implications. MostT Tauri stars are members of binary or multiple systems, and ina number of these it is believed that the disc and orbit planes aremisaligned. One particular system for which there is a resolvedimage of a disc which is misaligned from the binary plane is .M. Fragner & R.P. Nelson: Warped Accretion Discs in Binary Star Systems 25 HK Tau (Stapelfeldt et al. 1998). In this system, the disc radiusis estimated to be 105 AU, and the projected binary separationis approximately 3 times larger than this, suggesting that thedisc may be tidally truncated and undergoing strong interac-tion with the inclined companion star. The images do not showany signs that the disc is warped, which is consistent with our h = .
05 models, whose aspect ratio is probably appropriatefor T Tauri discs. Evidence for misaligned young binaries isalso provided by observations of precessing jets in star form-ing regions Terquem et al. (1999). In one particular example,Bally & Devine (1994) suggest that the jet which appears tobe launched by the young stellar object HH43* in the L1641molecular cloud in Orion precesses with a period of about 10 years. Applying models to a disc with outer edge of 50 AU sur-rounding a solar type star, the rigidly precessing disc modelspresented in this work precess with a period of between 3.5 –5 × years, depending on the inclination angle considered.Thus the observed precession is consistent with that driven bya binary companion, with parameters close to those adoptedin this work. The long alignment timescales that we find for h = .
05 discs suggests that a T Tauri disc will remain mis-aligned throughout its lifetime, such that jets launched fromthe disc inner regions will continue to precess for the durationof the time when jet launching is active.The thin h = .
01 disc models that we have computed aremost likely relevant to discs in X-ray binaries. The X-ray bi-nary systems SS433 Margon (1984) and Her X-1 Boynton et al.(1980) are two examples of systems thought to contain pre-cessing discs, with the evidence provided by the precessingjet of SS433 being particularly compelling. It is likely thatthe transfer of matter between the donor star and compact ob-ject in these systems will cause disc material to lie in the bi-nary orbit plane initially, and a tilting mechanism is requiredto move the disc out of this plane, such as the Bardeen-Petterson e ff ect (Bardeen & Petterson 1975), a misaligneddipole magnetic field associated with the central neutron star(Terquem & Papaloizou 2000), a disc wind (Schandl & Meyer1994) or radiation pressure (Iping & Petterson 1990; Pringle1996). If our results on the rapid realignment of highly twisteddiscs are correct, then this implies that a strong tilting mech-anism which operates on short timescales will be required tomaintain a misaligned disc which can be observed to precess.Finally, we briefly consider the implications of our resultsfor the formation of planets in misaligned binary systems. Theearly stages of planet formation are believed to involve thegrowth of planetesimals via mutual collisions. As the plan-etesimals grow in size, their interaction with the disc throughgas drag forces will decrease. Planetesimals orbiting at di ff er-ent radii will have their own precession frequencies, likely tobe di ff erent to that of the disc. Consequently, the growth ofthe planetesimals will eventually cause them to develop orbitswhich are increasingly inclined to the disc midplane. The e ff ectof this is likely to be an increase in relative motion betweenplanetesimals, which will a ff ect the outcomes of mutual colli-sions, and the rate at which planetary growth proceeds. A de-tailed examination of these issues will be presented in a forth-coming paper. References
Artymowicz P., Lubow S.H., 1994, ApJ, 421, 651Bally J., Debine D., 1994, ApJ, 428, L65Bardeen J.M., Petterson J.A., 1975, ApJ, 195, L65Boynton, P. E., Crosa, L. M., & Deeter, J. E. 1980, ApJ, 237,169Edwards S., Cabrit S., Strom S.E., Heyer I., Strom K.M.,Anerson E., 1987, ApJ, 321, 473Iping, R. C., & Petterson, J. A. 1990, A&A, 239, 221Katz J.I. 1980, ApJ, 236, L127Kley, W., & Nelson, R. P. 2008, A&A, 486, 617Kley W., Papaloziou J.C.B., Ogilvie G.I., 2008, A&A, 487, 671K¨onigl A., Ruden S.P., 1993,in Levy E.H., Lunine J., eds,Protostars and Planets III Univ. Arizona Press, Tucson, p.641Larwood, J.D., Nelson, R.P., Papaloizou, J.C.B., Terquem, C.,1996, MNRAS, 282,597Leinert C., Zinnecker H., Weitzel N., Christou J.,Ridgway S.T., Jameson R.F., Haas M., Lenzen R., 1993,A&A, 278, 129Lubow, S. H., & Ogilvie, G. I. 2000, ApJ, 538, 326Margon, B. 1984, ARA&A, 22, 507Nelson, R.P. & Papaloizou, J.C.B., 1999, MNRAS 309,929Ogilvie, G. I. 2000, MNRAS, 317, 607Papaloizou, J.C.B. & Lin, D.N.C., 1995, ApJ, 438,841Papaloizou, J.C.B. & Pringle, J.E., 1983, MNRAS, 202,1181Papaloizou, J.C.B. & Terquem, C., 1995, MNRAS, 274,987Petterson, J.A. 1977, ApJ, 214,550Pringle, J. E. 1996, MNRAS, 281, 357Schandl, S., & Meyer, F. 1994, A&A, 289, 149Shakura, N.I. & Sunyaev, R.A., 1973 A&A, 24 337Stapelfeldt, K. R., Krist, J. E., Menard, F., Bouvier, J., Padgett,D. L., & Burrows, C. J. 1998, ApJ, 502, L65Terquem C., Bertout C., 1993, A& A, 274, 291Terquem C., Bertout C., 1996, MNRAS, 279, 415Terquem C., Eisl¨o ffff