A&A 511, A77 (2010)

Evolution of warped and twisted accretion discs in close binary systems

M. M. Fragner - R. P. Nelson

Astronomy Unit, Queen Mary, University of London, Mile End Road, London E1 4NS, UK

Received 7 August 2009 / Accepted 2 October 2009

Abstract
Context. There are numerous examples of accretion discs in binary systems where the disc midplane is believed to be inclined relative to the binary orbit plane. Examples include the X-ray binaries Her X-1 and SS433, and the young stellar binary HK Tau. Under suitable physical conditions, such a configuration is expected to induce warping and rigid-body precession of the disc.
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, viscosity parameter , disc outer radius R, and binary inclination angle . We also aim to examine the conditions that lead to an inclined disc being disrupted by strong differential precession.
Methods. We use a grid-based hydrodynamic code to perform 3D simulations. This code has a relatively low numerical viscosity compared with the SPH schemes that have been used previously to study inclined discs. This allows the influence of viscosity on the disc evolution to be tightly controlled.
Results. We find that for thick discs (h=0.05) with low , efficient warp communication in the discs allows them to precess as rigid bodies with very little warping or twisting. Such discs are observed to align with the binary orbit plane on the viscous evolution time. Thinner discs with higher viscosity, in which warp communication is less efficient, develop significant twists before achieving a state of rigid-body precession. Under the most extreme conditions we consider (h=0.01, and ), we find that discs can become broken or disrupted by strong differential precession. Discs that become highly twisted are observed to align with the binary orbit plane on timescales much shorter than 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 precess rigidly. We also find that as h is decreased substantially, discs may be disrupted by strong differential precession, but for disc thicknesses that are significantly less (h=0.01) than those found in previous studies (h=0.03).

Key words: accretion, accretion disks - hydrodynamics - protoplanetary disks - binaries: close - binaries: general - methods: numerical

1 Introduction

Accretion discs are found in a number of astrophysical environments. 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 binary systems, where the peak in the distribution of separations occurs at approximately 30 AU (Leinert et al. 1993). Studies of the tidal interaction between a circumstellar disc and a binary companion indicate that tidal truncation occurs at a ratio of binary separation to disc radius for an equal mass binary (Kley & Nelson 2008; Artymowicz & Lubow 1994; Larwood et al. 1996). Given that 30 AU is smaller than typical disc sizes observed in T Tauri systems (Edwards et al. 1987), protostellar discs in many binary systems are often going to experience strong tidal effects.

There are observational indications that the disc and binary orbital plane are not always coplanar. The binary system HK Tau has been imaged, and shows compelling evidence for a disc whose midplane is misaligned with the orbit of the binary companion (Stapelfeldt et al. 1998). Terquem & Bertout (1996,1993) have suggested that the reprocessing of radiation from the central star by a tilted and warped accretion disc could account for the high spectral index of some T Tauri stars. A number of precessing jets have been observed in star forming regions, and the jet precession may indicate an underlying precessing disc from which the jet is launched (Terquem et al. 1999; Königl & Ruden 1993). In the case of close X-ray binaries, the objects SS433 and Her X-1 show evidence of precessing jets and/or accretion discs (Katz 1980; Petterson 1977; Boynton et al. 1980; Margon 1984).

An early study of the equations describing the behaviour of a thin non-coplanar viscous disc was undertaken by Petterson (1977). The resulting equations, however, lacked the property of conserving global angular momentum. This was addressed in subsequent work by Papaloizou & Pringle (1983), where they showed that warps in discs for which evolve diffusively on a timescale shorter than the usual viscous accretion by a factor , where is the dimensionless viscosity coefficient (Shakura & Sunyaev 1973), and h is the disc aspect ratio H/r. The opposite regime has been investigated by Papaloizou & Lin (1995) using linear perturbation theory. They demonstrated that the governing equation for the disc tilt changes from being a diffusion equation to a wave equation in this physical regime, showing that warps propagate as bending waves with a speed corresponding to half the sound speed, . An analytic study of warped discs was undertaken by Ogilvie (2000) in which the mildly non linear evolution of disc warps was examined.

The tidal perturbation of an inviscid disc by a companion on an inclined circular orbit has been studied by Papaloizou & Terquem (1995) using linear theory. They investigate the disc response to tidal perturbations that have odd symmetry with respect to the disc midplane and consider both zero and non-zero perturbing frequencies. They show that the zero-frequency (secular) mode leads to rigid precession of the disc if the sound crossing time through the disc is smaller than the differential precession time. Furthermore, they estimate that the timescale for damping the inclination is of the same order as the accretion timescale of the disc. A study of the twisting and warping of viscous discs in binary systems was undertaken by Lubow & Ogilvie (2000) using linear theory, and it was shown that discs whose outer disc radii are smaller than the truncation radius maybe unstable to tilting.

Numerical simulations of tidally interacting discs which are not coplanar with the binary orbit have been performed by Larwood et al. (1996) using SPH simulations. They found that the disc precesses approximately as a rigid body as long as the disc aspect ratio is large enough. For smaller aspect ratios their results suggest that the disc develops a modest warp, while for the disc becomes disrupted as a consequence of differential precession. Their results suggest that a disc can split-up into two distinct parts that behave independently. These SPH simulations also showed that the inclination evolves on the viscous timescale, as expected from Papaloizou & Terquem (1995).

The goal of the present study is to examine in detail the structure of misaligned accretion discs in close binary systems as a function of the important physical parameters: h, , the outer disc radius R, and the inclination angle . The binary companion 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 compared with the SPH schemes used to perform previous studies. This allows us to tightly control the disc viscosity, and therefore examine in more detail its influence on the disc evolution. We consider a broad range of disc parameters in which warps propagate either as bending waves or diffusively . As well as examining the quasi-steady disc structures which arise because of tidal interaction with the inclined companion, we also examine the conditions under which a disc becomes disrupted due to strong differential precession. In particular we are interested in examining whether the following criterion holds: a disc will achieve a state of rigid-body precession if warp propagation across the disc (either through bending waves or diffusion) occurs on a timescale shorter than the differential precession time.

The plan of the paper is as follows. In Sect. 2 we present the basic equations of the problem. In Sect. 3 we describe the numerical methods used in the simulations, and in Sect. 4 we discuss the linear theory of bending waves and calibrate the simulation code against calculations based on linear theory. In Sect. 5 we present our results for the tilted, tidally interacting discs, and in Sect. 6 we discuss our results and present our conclusions.

2 Basic equations

We consider the evolution of a thin, gaseous, viscous disc in orbit around a central star, where the disc is perturbed by a companion star orbiting in a plane which is not coincident with the disc midplane. For a range of physical parameters, it is expected that the disc will precess rigidly around the angular momentum vector of the binary system. The equations of continuity and motion for a viscous fluid in a precessing reference frame may be written as:

where

is the convective derivative, is the density, the velocity and P the pressure. The precession frequency is given by , and the disc angular momentum vector is assumed to precess around a vector pointing in the direction of . The gravitational potential is . We adopt a locally isothermal equation of state:

where the isothermal sound speed is defined by . Here h is the aspect ratio of the disc, H/r, and is the local keplerian velocity. is the viscous force per unit mass:

where is the viscosity stress tensor:

We adopt the standard ``alpha'' model (Shakura & Sunyaev 1973) to specify the kinematic viscosity . We solve the above equations numerically, using a system of spherical coordinates  , , ).

2.1 Orbital configuration

We work in a reference frame in which the origin of the coordinate system is fixed on the central star, and the secondary star moves on a circular orbit about the origin with position vector . The gravitational potential, , at any position vector is therefore given by:

where G is the gravitational constant, and and are the primary and secondary masses, respectively. The first term is due to the primary star, and the remaining two terms are due to the companion star. The last indirect term accounts for the acceleration of the origin of the coordinate system, which coincides with the location of the primary star. Note that there are no contributions to the potential due to the disc since we do not include disc self-gravity in our calculations, and the stars do not experience a gravitational force due to the disc. Consequently the binary orbital parameters remain unchanged.

Hence the secondary star feels the acceleration due to the central star and the indirect component of the potential. In a non-precessing reference frame, centred on the central star, whose x-y plane coincides with the orbital plane of the binary, the position vector of the companion in a circular orbit may be written as:

where D is the constant binary separation and is the binary orbital angular velocity.

3 Numerical methods

The system of equations described in Sect. 2 is integrated using the grid-based hydrodynamics code NIRVANA (Ziegler & Yorke 1997), adapted to solve the equations in a precessing reference frame. This code uses operator-splitting, and the advection routine uses a second-order accurate monotonic transport algorithm (van Leer 1977).

3.1 Boundary conditions

Periodic boundary conditions were applied in the azimuthal direction. At all other boundaries outflow conditions were adopted using zero-gradient extrapolation. Velocities at the radial boundaries were set to ensure that zero viscous stress occurs there.

3.2 Units

We adopt a system of units such that the unit of mass is equal to that of the central star, , the unit of radius is equal to the radius of the inner boundary of our computational domain, and we set the gravitational constant G=1. The unit of time then becomes , where is the keplerian angular velocity evaluated at r=1. When discussing simulation results, however, we will express time in units of the orbital period at r=10, the nominal outer disc radius in the simulations. Thus we describe a time interval of as ``one orbit''. Inclination and precession angles are described in units of degrees.

3.3 Reference frames and initial conditions

The numerical domain extends radially from r=[1,R], and azimuthally from . The outer radius, R, may differ from simulation to simulation. The meridional domain covers the range , where again may vary between the simulations.

We perform simulations in two different reference frames. For models in which the disc is expected to develop rigid precession, we solve the equations in a frame which precesses around the angular momentum vector of the binary. Given a sensible choice of the precession rate, this approach ensures that the disc midplane always stays close to the equatorial plane of the computational domain where , and allows simulations to be performed with large inclinations of the binary orbit. If evolved in a non-precessing frame, the disc would precess around the binary angular momentum vector, causing it to eventually intersect with the upper and lower meridional boundaries of the computational domain if the binary inclination angle is large. Adopting a precessing reference frame allows simulations to be conducted with relatively small meridional domains, even for large binary inclinations, thus substantially reducing the computational expense involved.

For models in which substantial differential precession of the disc is expected to arise, adopting a frame with a single precession frequency does not solve the problem described above, since some parts of the disc inevitably intersect the meridional boundaries. As these models are of significant interest from the astrophysical point of view (they will tend to be thin discs with large viscosity, similar to those which occur in X-ray binaries for example), we have computed them in a non precessing frame, but have been forced to use large meridional domains and smaller binary inclination angles.

We denote coordinates in the precessing frame as , while coordinates in the non-precessing binary frame are denoted . The coordinates in the code frame are denoted with , where the code frame can be one of either the precessing frame or the binary frame. The transformation of any vector from the binary into the precessing frame is given by:

with rotation matrices:

Thus to transform a vector from the binary into the precessing frame, we first rotate around the x-axis by an inclination angle and then rotate the resulting vector around the z-axis by an precession angle . Here is the precession rate of the precessing frame and t is the time elapsed. This transformation is particularly useful when setting up the initial conditions for the simulations which are performed in the non-precessing binary frame, since we assume that the binary orbit lies in a plane which is coincident with the equatorial plane of the computational grid, with the disc being inclined by an angle .

3.3.1 Model set-up in precessing frame

In the precessing frame the code midplane (defined to reside where ) coincides with the disc midplane initially, so that ( ). The initial density profile is chosen to satisfy hydrostatic equilibrium in the -direction, which yields the result:

where is the radial power law exponent of the density profile. Note that this expression has the limit:

which is the familiar form of the density profile for thin discs. The velocity in the direction is chosen to balance the forces in the radial direction:

Consequently, the azimuthal velocity profile is slightly subkeplerian because of the radial pressure gradient. The velocities in the r and direction are chosen to be zero initially.

We set in Eq. (2) to point in the same direction as the binary angular momentum vector:

where the last equality is derived using the transformation given by Eq. (8). The precession rate may be estimated using the expression:

where is the disc surface density. Hence the precession frequency is determined by the ratio of the integrated external torque exerted by the secondary to the integrated angular momentum content of the disc (Papaloizou & Terquem 1995). For the velocity and density profiles given by Eqs. (12) and (10), respectively, this evaluates to:

where is the orbital angular velocity of the disc at radius R. The binary orbit given by Eq. (7), expressed in precessing frame coordinates, is given by:

Thus an observer moving with the disc sees an increased binary frequency due to the retrograde precession of the disc (i.e. is negative).

3.3.2 Model set-up in binary frame

When working in the non-precessing binary frame we set the precession frequency in Eq. (2). Since the equatorial plane of the computational domain coincides with the binary orbit plane ( ), the disc is inclined with respect to the equatorial plane of the computational grid by an inclination angle . We can use the transformation defined by Eq. (8) at t=0 to relate the polar angles of the precessing frame to the polar angles of the binary frame. One particularly useful relation is given by:

which can be used to determine the initial density profile given by Eq. (10). The initial velocity is given by:

Using the relationship between the polar angles one can find an expression for that is written entirely in terms of binary frame angles. The different components in the radial, azimuthal and meridional direction are then given by:

where is given by Eq. (12) and by Eq. (17). It can be shown that this set of profiles satisfies the Euler equations given by Eqs. (1) and (2) with , , and , since in a spherically symmetric potential there is no preferred direction for the disc rotation axis.

3.4 Definition of precession and inclination angles

In order to follow the evolution of the disc structure we calculate the total angular momentum vector, , for disc material contained in each spherical shell in the numerical domain:

where is the radial grid spacing and is the angular momentum density:

and are the code unit vectors in the meridional and azimuthal directions, respectively. Note that when working in the precessing frame a term involving the precessional velocities should be added to this expression. In practice we find that this is of negligible magnitude, because of the slow rate of precession, and therefore the term has been omitted. For disc material located in a given spherical shell we locate the inclination angle, , (equal to at time t=0) between the disc and binary angular momentum vectors through:

where represents the binary angular momentum vector. The precession angle, , measured in the binary plane can be defined through:

where is an arbitrary reference unit vector in the binary plane. We choose , so that the initial precession angle by this definition. When working in the precessing frame has to be expressed in precessing frame coordinates and becomes time dependent.

When discussing the results of the simulations we use the the following terminology: a warp occurs if the angle becomes a function of radius in the disc; a disc develops a twist when  becomes a function of radius; a disc is said to be broken if either  or change discontinuously at some radius such that the disc separates into two independently precessing parts; a disc is said to be disrupted if varies smoothly by more than 180 degrees across the disc radius because of differential precession.

3.5 Calculation of torque contributions

Later in this paper we will be interested in measuring the different contributions to the changes in the precession and inclination angles for disc material located at different radii in the computational domain. As a first step we calculate the rate of change of the angular momentum vector associated with disc material located at radius r, which involves integrating over the individual spherical shells in the computational domain.

Using the continuity and momentum Eqs. (1) and (2) to express the velocity and density changes in Eq. (20), and assuming that the density vanishes at the meridional boundaries, we find the result:

where a dot denotes a time derivative. We have that

is the change due to the divergence of a radial angular momentum flux. The change due to viscous interaction with neighboring shells is given by:

Thus the z-component of the angular momentum vector can only be changed by the component of the viscous stress tensor, while the x and y components can also be changed due to radial shear of vertical motion via . The last term in Eq. (23) which induces a change of the angular momentum vector of disc material contained in a spherical shell is due to the gravitational interaction with the secondary star:

where the angular momentum change due to torques is given by:

and the torque density is given by:

When working in the precessing frame there is an additional term causing angular momentum change due to Coriolis and centrifugal forces:

with

In the following we are only interested in the torque components expressed in precessing frame coordinates. Thus all the torque components have been calculated now for simulations performed in the precessing frame. However, when working in the binary frame, we perform the above calculations and then project each of the torque components onto precessing frame coordinate axes to ensure a uniform approach to monitoring the torque contributions. Since the latter are time dependent, the additional term when working in the binary frame is given by:

This term accounts for the angular momentum change measured in the precessing frame that arise because of the precession of the frame.

3.6 Calculation of precession and inclination rates

We can now relate the above torques to the rate of change of the precession and inclination angles. When working in the precessing 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:

Note that the vector components L_X, L_Y and L_Z in the above expression, and in the remaining expressions in this subsection, should have hat symbols because they are measured in the precessing frame. We neglect these, however, in order to simplify the notation. For sufficiently rigidly precessing discs the deviation from the mean precession frequency will be small, and the disc midplane will remain close to the equatorial plane of the computational domain. In this case we can use the approximation:

With this approximation, Eq. (30) gives to first order:

Thus for sufficiently rigidly precessing discs the total inclination is the sum of the inclination of the precessing frame and a small deviation term which is proportional to the y-component of the angular momentum vector, measured in the precessing frame. Likewise the precession angle is the sum of precession due to the frame and a small deviation term that is proportional to the x-component of the angular momentum vector. The advantage of this approximation is that the precession and inclination angles become linear functions of the torque components. This allows us to write:

where each of the rates are calculated according to Eq. (32):

Hence, in circumstances where the approximation expressed in Eq. (31) holds, we can measure separate contributions to the inclination and precession rate coming from the radial flux (Eq. (24)), viscous stress (Eq. (25)) and gravitational interaction with the secondary star (Eq. (26)). When discussing contributions to differential precession rates later in the paper, we will omit the constant mean precession rate of the precessing frame 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 well it is able to model the propagation of free bending waves, and compare the results with linear theory. The linear theory of warps in accretion discs was initially investigated by Papaloizou & Pringle (1983). Their analysis is valid in a regime in which the dimensionless viscosity parameter (Shakura & Sunyaev 1973) lies in the range . In this regime, globally warped discs were found to evolve diffusively, with a diffusion coefficient larger than that associated with mass flow through the disc by a factor of , assuming an isotropic viscosity. Papaloizou & Lin (1995) investigated the regime in which , and showed under that, for low frequency modes, the governing equation for the disc tilt changes from being a diffusion equation to a wave equation, such that warps are communicated through the disc via propagation of bending waves. In a keplerian disc these are expected to propagate with little dispersion at half the sound speed, .

Nelson & Papaloizou (1999) used smooth particle hydrodynamic simulations to examine the propagation of free bending waves, comparing their numerical results with linear calculations. Overall, they found that the SPH code managed to model the propagation of bending waves quite accurately as long as the warp amplitude remains small such that linear theory holds. One drawback associated with using SPH for that problem, however, was the fairly large numerical viscosity associated with the scheme which is difficult to control and specify ab initio.

In this section we compare our numerical results with a linear calculation, adopting a similar set-up to the one used by Nelson & Papaloizou (1999). The equations describing the evolution of linear bending waves have been derived by Nelson & Papaloizou (1999). We do not reproduce the full derivation here, but rather mention some of the salient points associated with the derivation.

4.1 An initial value problem

Small-amplitude warps can be considered to be linear perturbations of a disc with midplane initially coincident with the (, ) plane. Taking the dependence of the perturbations to be , with m=1 for global warps, the Euler equations written in cylindrical coordinates can be linearised (Nelson & Papaloizou 1999; Papaloizou & Lin 1995). An analytical solution to the linearised equations corresponds to the case of a rigid tilt. In this case the velocity and pressure perturbations are given by:

where is the Keplerian angular frequency and is a small constant inclination. For large scale warps, it is expected that the local inclination varies on a length scale significantly larger than the disc thickness . It is then reasonable to 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 introducing the new quantities and . The resulting expressions are given by Nelson & Papaloizou (1999) in the form:

where

Equations (36) are a set of coupled differential equations which describe the evolution of the inclination as a function of r and t. Thus they can be used to follow the evolution of a bending disturbance from initial data. Below we use the time-dependent solution calculated in this way to test NIRVANA's ability to propagate bending disturbances.

Assuming a time dependence of the perturbations , with the wave frequency obeying the relation for slowly varying warps, one can derive algebraic equations for the perturbed quantities. These show that warps induce horizontal motions and vertical shear in the disc through terms proportional to , (Papaloizou & Lin 1995). Physically, these horizontal motions are generated by pressure forces in the disc which arise from the vertical misalignment of neighboring regions of the disc midplane. They are responsible for driving the bending wave, which propagates with approximately half the sound speed.

The main effect of viscosity in the disc is to damp these horizontal motions, reducing the efficacy of bending wave propagation. As the viscosity increases and , the bending disturbance begins to evolve diffusively on a time scale that scales inversely with (Papaloizou & Pringle 1983), where this latter quantity is now the diffusion coefficient associated with warp propagation in this regime. If the viscosity is small enough that its effect on the unperturbed quantities can be neglected, then the set of Eqs. (36) can be extended to include the effect of a small viscosity by the replacement (Papaloizou & Lin 1995):

If the viscosity is large, however, then the complete viscous stress tensor should be included in the equations of motion, making the analysis rather complicated, since there would then be accretion velocities in the unperturbed state.

4.2 Setting up initial data

A number of simulations using NIRVANA have been performed, in which warps with different amplitudes in discs with different parameters were set up and allowed to evolve. We compare these solutions with those obtained using the linear theory expressed in Eqs. (36) using the same set of initial conditions. The linear calculations employed a one dimensional finite difference scheme for which we used 1000 equally spaced grid cells. Starting with a Gaussian initial inclination profile:

with to be centred in the middle of a disc with inner radius r=1 and outer radius r=9, the set of Eqs. (36) was integrated forward in time. We chose the viscosity in these test simulations to be since we are mainly interested in the regime of propagating bending waves, and because the linear solution presented here becomes invalid for large viscosities. The NIRVANA code was set up in an inertial frame using the expressions (10) and (18) to specify the initial density and velocity, in combination with Eq. (17), with now being the disc tilt which is a function of r. In these simulations we used  cell in radius, 5 cells per pressure scale height in the meridional direction (and different numbers of vertical scale heights, depending on the initial warp amplitude, as described below) and  cells in azimuth. We used the same boundary condition as for the global simulations, which were described in Sect. 3.1.

4.3 Results

An example result is shown in Fig. 1 which has a disc aspect ratio of h=0.03 and maximum initial inclination angle . In this particular simulation the meridional domain contains 15 pressure scale heights ( ) in total. Figure 1 shows the inclination as a function of radius at different stages of the evolution. The solid blue line represents the NIRVANA solution, while the red dashed line shows the linear calculation.

Moving from the top left to top right panel, we can observe how the initial pulse begins to broaden. At time t=0.96 we see that the pulse has separated into an in-going and out-going bending wave. By time t=1.46 the pulses have reached the boundaries of the domain and have fully separated. It is clear that NIRVANA is able to capture the evolution of these low amplitude bending waves with a high level of accuracy.

In Nelson & Papaloizou (1999) the numerical SPH and linear solutions showed more substantial disagreement at the position of the initial pulse (r=5) at the final stage of the simulation, with the separation of the two pulses not being captured as accurately. This was probably due to the larger numerical diffusion present in the SPH code, though it should be remembered that the SPH simulations were performed using only 20 000 particles. The more accurate solution obtained with NIRVANA, on the other hand, is a clear indication that the numerical diffusion is smaller.

We also performed simulations for other values of the disc thickness h and maximum initial inclination angle . The number of pressure scale heights used in the meridional direction was scaled with the maximum inclination angle to avoid the disc interacting with the meridional boundaries. The results are shown in Fig. 2, and are presented at a stage of the simulations when the bending waves have propagated all the way to the radial boundaries. This corresponds to the final panel of Fig. 1. Since the bending waves propagate with half the sound speed, the evolution occurs over a longer time in the thinner discs (h=0.01), and over a shorter time in the thicker disc models (h=0.05). As can be seen from Fig. 2, the agreement between the numerical solution and the linear calculation is good in all the models presented here. This is particularly the case for those runs with lower initial maximum inclinations in Fig. 2 (top left and right panels), but the results are still in very reasonable agreement for the higher inclination cases in Fig. 2 (lower left and right panels).

We mention here that for even higher inclination angles NIRVANA behaves more diffusively, and the agreement between the linear and non linear solutions gets worse. As found by Nelson & Papaloizou (1999), when the warps become non-linear they generate horizontal motions that are of the order of the sound speed, creating shocks because of the symmetric form of the initial bending disturbance. The result is an effectively higher viscosity operating in the disc, which is clearly not accounted for in the linear calculations.

5 Tilted discs in binary star systems

We now discuss our results for the misaligned, tidally interacting discs. We set up the disc models according to the procedure outlined in Sect. 3, and details of the model parameters can be found in Table 1. The models are characterized by the disc aspect ratio, h, viscosity, , and initial inclination with respect to the binary orbit plane, . The disc outer edge, R, was chosen such that the disc is already very close to its tidal truncation radius, which is about one third of the distance between secondary and primary star, D. In some models the initial radius of the outer disc edge was varied also.

In this study, we are interested in examining the evolution of discs that fulfil the condition for wave-like warp propagation, , as well as discs that support diffusive warp propagation, . For the models in the wave regime we set , and for the runs in the diffusive regime we set .

The perturber is evolved on a circular orbit at a distance of D=30, and its mass is increased linearly from zero up to over a time interval of 4 orbits, during which time its distance is kept constant. Models 1-5 were performed in the precessing frame. If we use Eq. (15) to estimate the precession frequency, then we obtain a value of [degrees/orbit]. We find, however, that a better fit to test simulations is given by [degrees/orbit], which is the value used in the simulations. The reason is that as the disc gets tidally truncated it tends to precess slower, as can be seen from Eq. (15).

The resolution was chosen such there are 5 cells per pressure scale height in the meridional direction for all models. Models 1 and 3 incorporated 15 scale heights in the meridional direction. Models 2, 4 and 5 used 22.5 scale heights (to allow for a greater degree of twisting in the higher viscosity models where warp communication is expected to be less efficient). In the radial and azimuthal directions we used and  cells, respectively, for each of these models.

Table 1:   Table of runs.

Models 6 and 7 were performed in the non-precessing binary frame. Since the disc midplane is inclined with respect to the equatorial plane of the computational grid in these simulations, higher resolutions are necessary to avoid numerical errors. This is particularly true in the azimuthal direction because tilting the disc causes a component of the disc vertical structure to point along this direction. For a disc whose midplane coincides with the equatorial plane of the computational grid, pairwise cancelation of fluxes ensures conservation of angular momentum to machine accuracy. However, if the disc midplane is inclined with respect to the equatorial plane of the computational grid, the accuracy is limited by the advection scheme, and higher resolution in the azimuthal direction becomes necessary to avoid unphysical evolution of precession and inclination angles. Hence we used 600 cells in azimuth. In order to be able to resolve spiral waves for these extremely thin discs, we used 1056 cells in radius. The disc outer edge was chosen to be a bit smaller than in the thick disc runs, since highly pronounced non linear effects at the outer edge were seen in lower resolution test simulations (described as Models 6a and 7a in Table 1). In effect a strongly perturbed and narrow outer rim of gas became partially detached from the main body of the disc and developed a very different inclination and precession angle evolution in Model 6a, an effect not observed in Model 6 with a smaller initial outer radius. This phenomenon is discussed in more detail later in the paper. In order to accommodate an inclined disc with an inclination of we used 60 pressure scale heights (17.03 degrees) in the meridional direction, with 5 cells per pressure scale height.

We now describe the results of these models, where we have ordered our discussion in terms of the disc aspect ratio, h.

5.1 Models 1 and 2

A column density plot for the disc in Model 1 is displayed in the upper panels of Fig. 3 corresponding to the end stage of this simulation. The left panel displays a projection of the disc column density onto the (, ) plane of the binary reference frame, and the right panel displays a projection onto the binary (, ) plane. At this stage the disc has precessed rigidly to  degrees, so the disc appears edge on in the (, ) plane and face on in the (, ) plane. The precession of 180 degrees is evident, as the disc angular momentum vector now has a negative component, whereas it had a positive component in its initial state, as can be seen from the transformation given by Eq. (8) at t=0. Another result of the interaction with the secondary star is the formation of spiral waves, which are apparent in Fig. 3.

Figure 4 shows the evolution of the precession () and inclination () angles, respectively, for disc material orbiting at different radii. These have been calculated using Eqs. (22) and (21) respectively, where has been calculated using Eq. (19), which we have integrated over radial shells of width . Thus the red dashed line in Fig. 4 corresponds to disc gas between r=1 and r=2, the blue solid line to disc gas between r=8 and r=9, and the black dotted line to disc gas between r=4 and r=5. The solid black line shows the precession and inclination angles of the entire disc. It is apparent that the lines are almost indistinguishable in Fig. 4 (upper left panel), showing that the different disc parts in Model 1 precess at the same uniform rate, such that the disc is in a state of rigid body precession with almost no twist. This behaviour is expected from linear theory (Papaloizou & Terquem 1995). In Table 1 we show the warp propagation timescale and the differential precession timescale, . In the bending wave propagation regime (), , since warps propagate with half the sound speed in this regime. In the diffusive regime () . The warp propagation timescale should be compared to the differential precession timescale . If it is expected that the disc will precess as a rigid body, while if the disc may become disrupted as a result of differential precession. Hence, we see that the rigid precession of the disc in Model 1 agrees with expectations. Our Model 1 is very similar to Model 1 in Larwood et al. (1996), which also shows rigid body precession. The inferred precession period in our Model 1 is about 146 orbits, which is a little larger than the period given by Eq. (15). The difference arises because the disc is tidally truncated and shrinks slightly, so its precession rate decreases.

Looking at the upper right panel of Fig. 4, we observe that the inclination for Model 1 changes only slightly during the simulation. The oscillation seen in the inclination rates is driven by the secondary star. Since the orbital plane is inclined with respect to the disc midplane, the vertical component of the gravitational force due to the secondary star causes the disc to perform a forced oscillation with a frequency equal to twice the binary frequency. Papaloizou & Terquem (1995) argued that the disc is expected to align with the binary orbit on a time scale similar to that over which binary torques cause the total angular momentum of the disc to evolve. If the viscously-induced outward expansion of the disc is balanced by tides due to the companion, then disc-alignmnent should occur on the viscous evolution time, . For Model 1 we estimate orbits. Extrapolating the disc evolution shown in Fig. 4 gives an alignment time of 2025 orbits, in good agreement with the viscous time scale.

Lubow & Ogilvie (2000) present images which represent disc shapes for a variety of disc parameters, obtained using a linear analysis of tilted discs in binary systems. They introduce a parameter , which is a measure of the disc thickness, and is roughly equivalent to . Their model a displayed in their Fig. 7 has , and D/R=0.3, similar to our model 1. Although a detailed comparison is not possible, it appears that there is general agreement between linear theory and our non-linear simulation, since both show almost no distortion of the disc due to twisting or warping.

Column density plots for Model 2 are displayed in the lower panels of Fig. 3. Warp propagation in this higher viscosity model ( ) is expected to be diffusive, and therefore less efficient than for Model 1 (as shown by the warp propagation times listed in Table 1). This run was also evolved until the precession angle reached -180 degrees. The high viscosity operating in the disc leads to the damping of the spiral waves before they can propagate very far, and thus they are not apparent in the images.

The precession angles for Model 2 are displayed in Fig. 4 (lower left panel). We observe that this more viscous disc develops a small differential twist before it attains a state of rigid precession. This suggests that the disc needs to develop a slightly more distorted shape in order for internal stresses to counterbalance the companion-induced differential precession when warp communication is less efficient. Nonetheless, rigid body precession is expected in this case, and the numerical results are in agreement with this expectation.

Examining the inclination evolution for Model 2 shown in the lower right panel of Fig. 4, we see it is very similar to that observed for Model 1. The viscous time scale in Model 2 should be approximately 4 times shorter than for Model 1, but there is no indication that the inclination evolution is faster in this case. We note that the inclination evolution rates do not appear to have reached final steady values for either Models 1 or 2 by the end of the simulations, such that an accurate assessment of the longer term evolution rate is not possible. This may be because the outer disc radius has not had time to equilibrate during the simulations as they have only run for a fraction of the global viscous timescale of the disc. We are able to conclude, however, that the inclination evolution occurs over timescales much longer than the precession time, and appears to be consistent with evolution on the global viscous time of the disc.

5.2 Models 3, 4 and 5

As shown in Table 1, Models 3, 4 and 5 all have h=0.03. Model 3 has , and so warp propagation occurs via bending waves, whereas Models 4 and 5 have such that warps propagate diffusively. Table 1 also shows that the warp propagation for Models 4 and 5 is expected to be considerably less efficient than for Model 3. Models 3 and 4 had inclinations of , whereas Model 5 had an inclination of , which induces a more rapid precession of the disc, and therefore has a greater tendency to twist the disc up.

At the outset of this project it was our intention to run all models until they had precessed by . For models 3 and 4, however, we found that the discs become eccentric on a timescale of approximately 60 orbits (16 binary orbits), leading to undesirable changes in the disc structure due to our use of an open inner boundary condition. Thus we did not evolve these models until they had precessed by . We note that the timescale for the observed disc eccentricity growth corresponds closely to that obtained by Kley et al. (2008). Model 5 did not suffer from this problem because the lower inclination angle induced more rapid precession, such that this model could be evolved until it had precessed by , but prior to growth of significant eccentricity. The study of disc eccentricity goes beyond the scope of this paper, but the results for models 3 and 4 suggest that very long term integrations may lead to the formation of inclined and eccentric discs. The fact that we did not evolve models 3 and 4 through of precession means that the upper and middle panels of Fig. 5 display images which correspond to of precession.

It is clear from these panels that the disc in Model 3 has precessed as a rigid body, with almost no twist apparent. Consequently, after precessing through the disc appears edge on when projected into the (, ) plane and face on when projected into the (, ) plane. The middle right panel of Fig. 5, however, shows that the disc in Model 4 has developed a small twist due to the inner disc not having precessed by at the same moment when the outer disc has. The lower panels of Fig. 5 show column density plots for Model 5, which was evolved until the outer disc had precessed by (assisted by the fact that the precession is faster in this case). It is apparent from these panels that the disc has developed a small twist since the inner disc has not precessed a full by the end of the simulation.

The precession angles for Model 3 are plotted in the top left panel of Fig. 6. The blue solid line corresponds to disc material orbiting between radii r= 8-9, the red dashed line radii between r= 1-2, and the black dotted line to radii between r= 4-5. A black solid line is also plotted showing the precession angle integrated over the whole disc. The fact that these lines effectively lie on top of each other shows that the disc precesses as a rigid body with essentially no twist. Rigid-body precession is expected in this case since the warp propagation time is less than the differential precession time (see Table 1).

The inclination angles for Model 3 are plotted in the top right panel of Fig. 6. The line styles and colours correspond to the same disc annuli as described above. It is clear that the inclination evolution in this case is very slow. The global viscous evolution time for this disc is approximately 10 000 orbits, and a linear extrapolation of the inclination evolution shown in Fig. 6 gives a time of 13 500 orbits for the disc to fully align with the binary orbit. Clearly the disc inclination is evolving on the viscous time in this case.

The evolution of the precession angles for Model 4 are shown in the left middle panel of Fig. 6. The line styles and colours correspond to disc annuli as described above. We see that this disc develops a well-defined twist which is set up within the first 5-10 orbits due to the inner and outer parts of the disc undergoing differential precession, with the outer disc precessing faster than the inner disc. The magnitude of the twist is approximately . Thereafter the disc is seen to precess as a rigid body, maintaining the same twisted shape. We note that this disc is expected to achieve rigid body precession from consideration of the warp propagation and differential precession times (see Table 1).

The inclination angles for Model 4 are shown in the right middle panel of Fig. 6. We see that the disc in this case has developed a small warp, with the inner and outer discs having a difference of inclination of about . We also observe that the inclination is evolving more rapidly than is the case for Model 3. The global viscous evolution time for this disc is approximately 1500 orbits, and extrapolating the inclination evolution observed in Fig. 6 gives an estimated time for alignment of 1000 orbits, close to the expected value.

It is interesting at this point to compare our results for Models 3 and 4 with Model 9 presented by Larwood et al. (1996), which had comparable parameters (h=0.03, D/R=3, 45 degree inclination). Although it is difficult to know the precise value of the viscosity operating in the SPH simulations, it is likely that the effective value lies between the values used in our Model 3 ( ) and Model 4 ( ). Nonetheless, we find qualitatively different behaviour, with our models quickly achieving a state of rigid body precession, with a small twist and warp which vary smoothly across the disc. Model 9 of Larwood et al. (1996) shows significant disruption of the disc, which effective breaks discontinuously into two distincts parts which become mutually inclined by approximately . The origin of this discrepancy is unclear, but one possibility is that SPH simulations using 20 000 particles are only marginally able to resolve the vertical structure of a disc with h=0.03. This may have the effect of reducing the efficacy of warp propagation below that which is expected from consideration of the disc physical parameters.

We now turn to Model 5. The precession angles for this run are displayed in the lower left panel of Fig. 6. Given that the precession frequency, given by Eq. (15), is proportional to , the disc precesses faster in this simulation than in Models 1-4, and this faster precession can be observed in Fig. 6. The lower left panel of this figure shows the precession angles as a function of time for disc annuli located between r= 1-2, r= 4-5, and r= 8-9, with the linestyles and colours being the same as described above. The larger precession frequency induces greater differential precession in the disc, and the consequence of this is that Model 5 shows a slightly larger twist ( ) than Model 4 ( ). Also plotted in this panel are precession angles that would be observed if the 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):

As can be seen from the figure, the disc parts tend initially towards their free particle rates, such that the outer disc precesses faster than the inner disc. Information about this differential precession is communicated across the disc, and internal stresses are established which cause the inner disc precession to speed up and the outer disc precession to slow down, such that precession at a uniform rate occurs after approximately 10 orbits. The warp communication and differential precession times given in Table 1 for this model predict rigid-body precession, as observed. The emergence of internal stresses as the disc twists up is illustrated in Fig. 7, where the different physical contributions to the precession rate are displayed for Model 5. These have been calculated according to the procedure outlined in Sect. 3.6. Here the different linestyles correspond to the same disc annuli described above (solid blue line: r= 8-9; dotted black line: r= 4-5; dashed red line: r= 1-2). From the lower left panel of Fig. 7, we can observe that the gravitational interaction with the secondary causes a negative precession rate for the outer annulus, and a positive precession rate for the inner annulus, and these rates correspond to the free particle precession rates. Since the frame in which these rates were measured precesses in a retrograde sense with a mean frequency , this means that the inner annulus would lag behind while the outer annulus would lead the precession of the reference frame, if the companion's gravity was the only contributor to the total precession rate. In fact, this is what is observed until about t=10 orbits, during which time the disc precesses differentially (see Fig. 6, lower left panel). As the disc starts to develop a twist due to differential precession, misalignment of disc midplanes as a function of radius causes pressure forces to generate horizontal shear motions that drive the propagation of bending disturbances, and these communicate information about precession angle changes across the disc. The resulting internal stresses are clearly proportional to the degree of twisting, allowing the disc to achieve a state of rigid body precession once a sufficiently twisted disc has been established. At this point the disc structure becomes quasi-static in the precessing frame as each disc annulus precesses at the same rate. This sequence of events is demonstrated by the upper left panel of Fig. 7, which shows 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 almost completely compensates the effect of gravity, such that the total precession rate (relative to the precession of the reference frame) is zero for all disc annuli, as can be seen in the lower right panel of Fig. 7, which shows the sum of the contributions to the precession rate. This demonstrates that the disc manages to redistribute the angular momentum such that rigid body precession is achieved. Once again we note that in this figure we have subtracted the mean precession rate of the precessing frame.

The effect of viscosity as calculated in Sect. 3.6 is depicted in the upper right panel of Fig. 7. We observe that the effect of viscous friction between differentially precessing radially adjacent disc shells causes the disc to precess more rigidly. However, the dominant effect of the viscosity is the damping of the horizontal shear motions induced by the twist. Since these shear motions are responsible for driving the bending disturbances, the viscosity will lead to weakened communication, and the redistribution of angular momentum due to radial advective fluxes shown in the upper left panel of Fig. 7 will be suppressed. This effect is dominant and opposite to the one seen in the upper right panel of Fig. 7.

The inclination angles for Model 5 are shown in the lower right panel of Fig. 6. It is clear that the disc develops a small warp, with the inner disc having higher inclination than the outer disc. The rate of global inclination change is similar to that observed in Model 4, suggesting alignment of the disc in the global viscous evolution timescale.

We plot the different physical contributions to the rate of inclination change for Model 5 in Fig. 8. In Fig. 9 we also display the time averaged contributions to the rate of inclination change as a function of radius, where the time averaging was performed over the full length of the simulation. As may be observed in the lower left panel of Fig. 8, and even more clearly in Fig. 9 (blue dotted line), the companion's gravity induces a globally negative net inclination change, leading to coplanarity of the entire disc on the viscous time scale, in agreement with the lower right panel of Fig. 6. This is consistent with the findings of Papaloizou & Terquem (1995) that the zero-frequency (secular) gravitational interaction will lead to coplanarity. We note that because of the time averaging, this is the dominant gravitational contribution that we pick up in Fig. 9 (blue dotted line).

More interestingly, we can observe from the upper left panel of Fig. 8, and in Fig. 9 (red dashed-dotted line), that the contribution from the radial fluxes not only counteracts the effect of gravity, it actually overshoots slightly, causing the inner disc to become more inclined than the outer disc. This change of inclination is seen in the upper right panel of Fig. 6. Thus the communication of bending disturbances on the one hand leads to rigid precession of the disc, and on the other hand causes the disc to become slightly warped. Viscosity tends to reduce the amount of warping, such that it acts to flatten the differential inclination profile across the disc radius, as can be seen from the green dashed line in Fig. 9. Thus the total average inclination rate in Fig. 9 (dashed-triple-dotted line) is negative for the entire disc, while the magnitude is slightly bigger for the outer disc, with the consequence that the disc has developed a slight warp and tends to become coplanar with the orbital plane on the viscous timescale, as seen in the lower right panel of Fig. 6.

5.3 Models 6 and 7

We now discuss Models 6 and 7, for which the discs have aspect ratio h=0.01, an inclination angle of , and were performed in a non precessing frame. Model 6 has , so warps propagate via bending waves in this case, and Model 7 had , so warps propagate diffusively in this model. We see from Table 1 that differential precession is expected to disrupt the discs in Model 7, as the warp propagation time is longer than the differential precession time. Rigid-body precession is expected for Model 6.

Column density plots showing the disc in Model 6 are presented in Fig. 10 for different times during the evolution. The top panels show the initial state of the disc. The disc appears edge on in the ( , ) plane and face on in the ( , ) plane. As we move down to the next panels the outer edge of the disc has precessed by 45 degrees at time t=10.4, as can be seen in Fig. 11 (upper left panel). The inner region of the disc, however, has only precessed by about 12 degrees at this stage and so appears almost edge-on in the left panel. In the third panels down the outer disc has precessed by 90 degrees at t=25.7 and appears edge on in the ( , ) plane and face on in the ( , ) plane. The inner disc has precessed by approximately 70 degrees by this time, and so does not appear edge-on in either projections. The bottom panels show the disc structure at the final output of the simulation which was halted at t=44 orbits, by which time the outer disc had precessed by approximately .

The precession angles for Model 6 are shown in the upper left panel of Fig. 11. The blue solid line corresponds to disc material orbiting in an annulus between r= 7-8, the black dotted line corresponds to the annulus between r= 4-5, and the red dashed line corresponds to the annulus between r= 1-2. As with Models 4 and 5, we see that the disc inner and outer annuli begin to precess at the local free-particle rates, and the disc develops a significant twist which reaches a maximum amplitude of approximately . As the disc becomes increasingly twisted, however, internal stresses are established which cause the inner disc precession rate to increase and the outer disc precession rate to decrease. Eventually the disc reaches a state of rigid body precession after a time of about 10 orbits. After this time we see that there is a long term readjustment of the degree of twist such that by t=44 orbits the twist angle has reduced from to .

In Fig. 12 we display the time integrated contributions to the precession rate as a function of radius for Model 6, following the procedure described in Sect. 3.6. Note that the average (negative) precession rate of the disc has been subtracted in this figure, such that a positive value corresponds to precession that will lag that of the overall disc, and a negative value corresponds to precession that leads. The contribution from the companion's gravity is shown by the blue dotted line, and this is very similar to the free-particle precession rate which is plotted using the black solid line. The red dashed-dotted line shows the contribution from the pressure-induced radial flux of angular momentum, and this is seen to almost exactly counterbalance the effect of gravity, showing that it is the pressure-induced radial motion in the disc which is responsible for halting the differential precession and inducing rigid body precession. The contribution from the disc viscosity is shown by the green dashed line, and is seen to be negligible.

The inclination angles for Model 6 are shown in the top right panel of Fig. 11, where the line styles and colours correspond to the disc annuli described above. As with Model 3, we see that the disc develops a warp such that the inner disc has a larger inclination to the binary plane than the outer disc. Indeed, the inclination of the inner disc is seen to increase initially, before decreasing after approximately 20 orbits as the disc as a whole aligns with the binary orbit plane. As discussed previously, disc alignment is expected to occur on the viscous time scale, which for Model 6 is approximately orbits. It is obvious from Fig. 11, however, that the disc in our simulation is not aligning on such a long time scale, but will actually align within 650 orbits. Possible reasons for this enhanced alignment rate are discussed below.

Figure 13 shows the time integrated contributions to the rate of inclination evolution as a function of radius for Model 6. We see that the companion's gravity causes the inclination to decrease in the disc inner and outer regions, and in the inner regions the radial flux of angular momentum opposes this and increases the inclination during the simulation. Indeed the radial flux contribution causes the inclination angle to grow there, as observed in the lower right panel of Fig. 11. The outer parts of the disc experience negative contributions from both the companion's gravity and the pressure-induced radial flux, driving alignment of the disc on the relatively short timescale described above. It is clear that the disc viscosity plays a negligible role in driving the inclination evolution (green dashed line).

We may compare our model 6 with model c in Fig. 7 of Lubow & Ogilvie (2000), which had , and D/R=0.3. Detailed comparison between the results presented in Lubow & Ogilvie (2000) and our model 6 is not possible, but we note that there is general agreement in the prediction that a disc with and will show significant distortion due to twisting and warping. A detailed comparison between our non linear simulation results and the predictions of linear theory will be presented in a future publication.

We now discuss the results from Model 7. In Fig. 14 we display column density plots for this model. The top panels show the initial state. The second panels down show the disc at time t=20.5. At this stage the outer disc edge has precessed by 90 degrees and appears edge on in the ( , ) plane, and face on in the ( , ) plane. By contrast, the inner disc has not precessed very much and still appears almost edge on in the ( , ) plane and face on in the ( , ) plane. At time t=44 orbits the outer disc has precessed by 180 degrees and appears edge on in the ( , ) plane and face on in the ( , ) plane. The inner disc, however, has only precessed by about 70 degrees by this time, and so does not appear edge on in either projection. The bottom panels show the final stage of the simulation when the outer disc has precessed by about 220 degrees. It is clear from these images that the disc in Model 7 has adopted a highly twisted shape, but one in which the precession angle for different disc annuli varies very smoothly across the disc.

The precession angles for Model 7 are displayed in the lower left panel of Fig. 11. We observe that the disc is in a state of differential precession at the end of the simulation, and the twist between inner and outer disc annuli is approximately 110 degrees at this stage. Although this twist is large, the actual rate of differential precession at the end of the simulation is very small, as the precession rates of the inner and outer disc have converged during the run. Indeed, extrapolation of the precession angles in Fig. 11 indicates that rigid body precession will be achieved at a time equal to t=56.6 orbits, at which point the twist amplitude will be . This suggests that even though the warp propagation time is much smaller than the differential precession time for Model 7 (see Table 1), internal stresses can be set up within the disc that cause uniform precession once the disc becomes highly twisted. Adopting a definition of disc disruption which requires the twist to become greater than 180 degrees, we suggest that the disc in Model 7 will not be disrupted, but will instead adopt a highly twisted but uniformly precessing configuration.

The inclination angles for Model 7 are shown in the lower right panel of Fig. 11, where the lines styles correspond to the disc annuli described above. We can see that the disc in this case develops a small warp with the difference between the inclinations of the inner and outer disc being . We also note that it appears that the integrated tilt of the whole disc is actually smaller than for any of the individual disc annuli, but this is an effect of averging over a significantly twisted disc. As observed for Model 6, we see that the inclination of the disc evolves rather quickly. The global viscous evolution time for Model 7 is 15 000 orbits, whereas extrapolation of the inclination angles for either the inner or outer disc suggests that they will completely align after 100 orbits.

What is the explanation for this rapid alignment observed for Models 6 and 7? One possibility that we have explored is that numerical diffusion may cause a tilted disc to align with the equatorial plane of the computational grid, since the advection routine does not conserve total angular momentum when the disc azimuthal velocity is not directed along the azimuthal direction of the computational grid. Lower resolution test calculations performed for a tilted disc without a companion indicate that this effect can cause eventual disc alignment, but at a rate which is at most 7 times smaller than observed in Models 6 and 7. This suggests that numerical diffusion is not the cause of the rapid alignment.

Another possibilty arises when we consider the evolution of a rigidly precessing, tilted disc which is in a precessing reference frame whose precession frequency is equal to that of the disc. Consider first the situation in which the angular momentum vectors for all disc annuli lie in the y-z plane of the binary frame initially (i.e. the z direction of the precessing frame lies initially in the y-z plane of the binary frame, pointing along the positive y and z directions). A torque exerted on disc annuli in the y direction, (a ``y-torque''), will cause the disc inclination to change. A torque exerted in the x direction (an ``x-torque''), will cause precession at a rate which differs from the precession rate of the frame. A positive x-torque causes precession to occur at a rate which is faster than the frame, and a negative x-torque causes slower precession. This general picture is expressed mathematically by Eqs. (34), and is shown diagramatically in the left hand image of Fig. 15. For a disc which precesses uniformly without any twist, then the companion's gravity will exert positive x-torques on the outer disc annuli, and negative x-torques on the inner disc annuli, as it tries to cause the disc to precess differentially. Internal stresses, however, will exactly counterbalance these x-torques to maintain uniform precession. In this simple scenario there is no torque exerted on the disc in the y direction to modify its inclination.

Consider now the case of a highly twisted disc which is being forced to precess uniformly, and whose annuli all have the same inclination. Gravity will again try to induce differential precession, but now the torque exerted on disc annuli will have both an x and a y component due to the twist. Resolving the vectors associated with the internal stresses required to enforce uniform precession shows that their combined x-torques operate in opposite directions such that they cancel when integrated over the disc. The combined y-torques, however, operate in the same direction, implying that the internal torques generate a net y-torque on the disc. Conservation of angular momentum obviously prohibits this from arising, suggesting that a y-torque must operate, whose effect 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 precession timescale in order to maintain uniform precession, this argument suggests that a highly twisted disc may align on the precession time, as is observed for Models 6 and 7. Interestingly, the resultant y-torque that arises in this picture is proportional to , where is the twist amplitude, and the rate of inclination evolution observed in Models 6 and 7 is in agreement with this scaling.

If the above argument is correct, then it suggests that discs with modest levels of twist will align on the global viscous timescale of the disc, but highly twisted discs may align on the much shorter precession timescale. This obviously has significant implications for astrophysical systems which may be expected to harbour highly twisted discs, such as the X-ray binaries Her X-1 and SS433.

5.3.1 Models 6a and 7a - broken or disrupted discs?

We now briefly discuss the results of two lower resolution models which were run during an early stage of this project: Models 6a and 7a. These were identical to Models 6 and 7 except that the number of grid cells in the radial direction was , and the disc outer radii were located at R=10 instead of R=8. Extending the disc radius outward has the potential effect of increasing the rate of differential precession across the disc, provided that tidal truncation does not cause the disc to shrink back down to R=8.

Figure 17: Evolution of precession ( left panels) and inclination angles ( right panels) for Models 6a and 7a. The blue solid lines correspond to disc annuli of unit radius lying between r=7 and r=10. The solid blue line denotes the annulus lying between r=9 and r=10. The red lines correspond to disc annuli of unit width lying between r=1 and =7. The solid red line denotes the annulus lying between r=1 and r=2.

Open with DEXTER

Column density plots for Model 6a are shown in the top two panels of Fig. 16, corresponding to an evolution time of t=22.8 orbits. We can see that there has been some tidal truncation of this disc, and the outer disc edge is significantly distorted in the images. This arises because a narrow outer rim, running between radii and , has effectively detached itself from the main body of the disc, and has evolved to have significantly different precession and inclination angles compared to the main body of the disc.

The inclination and precession angles for Model 6a are displayed in the upper left and right panels of Fig. 17, respectively. The precession angle of the annulus which lies between radii r= 1-2 is shown by the thick red solid line, and the remaining red dotted lines show the precession angles for disc annuli of unit width out to r=7. The solid blue line corresponds to the disc annulus lying between r= 9-10, and there are two additional blue dotted lines for annuli lying between r= 7-8 and r= 8-9. It is clear that the disc has separated into two parts discontinously at a radius between r=7 and r=8, with the detached outer rim having precessed significantly faster than the inner disc. The inclination angles plotted in the upper right panel show that the disc becomes significantly warped, with the outer rim tending to align with the binary and the inner disc actually increasing its inclination. Although the change in inclination appears to occur fairly smoothly across the disc, it should be noted that the disc has developed a fairly large warp of .

The behaviour of Model 6a is reminiscent of the broken disc presented by Larwood et al. (1996) (their broken Model 9 had h=0.03, inclination of , D/R=3) except that the discontinuous break in Model 6a occurs near to the edge of the disc. It would appear that Model 6a breaks because the larger disc radius induces a larger rate of differential precession, and also possibly because warp propagation is retarded in the outer disc regions which are subject to non linear density structures and shocks because of strong interaction with the binary companion.

Column density plots for Model 7a are shown in the lower two panels of Fig. 16, corresponding to a run time of t=27 orbits. Here we can see that the disc has become very highly twisted (the twist amplitude in the figure is ), but does not show any of the discontinuous behaviour observed in Model 6a. Instead the twist varies very smoothly across the disc. The precession and inclination angles for Model 7a are shown in the lower left and right panels of Fig. 17, respectively. The solid red line shows the precession angle of the annulus located between r= 1-2, and the solid blue line corresponds to the annulus between r= 9-10. By the end of the simulation the twist amplitude is , and the precession rates of inner and outer disc remain very different such that this disc will differentially precess through . So, by our working definition of disc disruption, this disc will be disrupted because of differential precession. The inclination angles show that the disc also becomes significantly warped, with the outer disc aligning with the binary orbit more rapidly than the inner disc.

To summarise: Model 6a and 7a both show examples of discs which break or are disrupted because of differential precession. In the case where warps propagate via bending waves, the disc breaks discontinuously into two independently precessing parts, but which themselves remain as coherent structures which precess as rigid bodies. In the case where warps propagate diffusively, viscosity is able to smooth out any discontinuous behaviour, and instead the disc becomes disrupted because the twist exceeds 180 degrees. In all probability the disc annuli in Model 7a will continue to precess differentially in the longer term, substantially destroying any coherent disc-like structure. The larger disc radius in this case, when compared with Model 7, causes the differential precession rate across the disc to be larger, and thus explains why Model 7a shows disruption, whereas Model 7 results in a highly twisted disc which is able to attain a state of rigid-body precession.

6 Conclusions

We have presented non linear simulations of accretion discs in close binary systems in which the binary midplane is misaligned from the binary orbital plane. Previous work on this problem, which employed linear theory, semi-analytic techniques and SPH simulations (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Larwood et al. 1996), suggests that under suitable conditions the disc should become mildly warped and precess as a rigid-body around the angular momentum vector of the binary system. The main aim of the present study is to examine in detail the structure and evolution of discs in misaligned binary systems, subject to variation of the important physical parameters. We used the grid-based code NIRVANA to perform the numerical simulations, and in principle this should have the advantage over previous SPH simulation studies that the disc viscosity is a well defined and controllable parameter.

Following Nelson & Papaloizou (1999), the propagation of linear bending waves was used as a means to calibrate the code and test its ability to propagate warps. Direct comparison with linear calculations show that the code can propagate bending waves with a high degree of accuracy for a range of disc models and warp perturbation amplitudes.

In our study of discs in misaligned binary systems, a number of different models were considered. Discs with aspect ratios h=0.05, 0.03 and 0.01 were studied, and for each of these values models were computed with viscosity parameters and . The smaller value of allows bending waves to propagate, whereas the larger value causes warps to evolve diffusively.

It is expected that discs for which the warp propagation time, , is shorter than the differential precession time, , will undergo rigid-body precession. Discs which do not satisfy this criterion may be disrupted by strong differential precession. All our models with h=0.05 and 0.03 have , and the simulations show that each of these discs attains a state of rigid-body precession. Discs with lower viscosity maintained a precessing structure with almost no discernible twist and warp, due to the highly efficient warp propagation associated with bending waves. Discs with larger viscosity propagate warps less efficiently, and these models underwent a short period of differential precession, setting up a twist in the disc, prior to attaining a state of rigid-body precession. Our analysis shows that internal stresses are established in the disc as it becomes twisted, which transport angular momentum across its radius and cause it to precess rigidly. Models in which warp propagation is the least efficient develop the largest twist, since the internal stresses which hold the disc together against differential precession are proportional to the twist amplitude. Examination of the simulations shows that the dominant contributor to these internal stresses are pressure-induced radial angular momentum fluxes which are driven by local misalignment of disc midplanes. When , the final amplitude of the twist was just a few degrees for the h=0.05 disc, and for the h=0.03 disc it was between 12-15 degrees, depending on the binary inclination. The low viscosity models showed essentially no twist.

In addition to becoming twisted, a disc may also become warped. We find that for all models with h=0.05 and h=0.03the degree of warping is very modest, with differences in inclination across all disc annuli being less than 1 degree. It is expected that these discs will also align with the binary orbital plane on the global viscous evolution time of the disc (Papaloizou & Terquem 1995; Larwood et al. 1996), and all simulations with h=0.05 and 0.03 are fully consistent with this expectation.

We also considered a number of thin disc models with h=0.01, where we not only varied the viscosity but also the outer radius (taking values of either R=8 or R=10). The low viscosity model with R=8 developed a significant twist of between 20- before achieving a state of solid body precession, as expected since for this model. A model run with and R=8 was predicted to be disrupted due to differential precession since in this case. During the simulation it developed a very large twist , but as the disc become increasingly twisted the rate of differential precession decreased. At the end of the run the disc was experiencing very slow differential precession, and extrapolation of the rate of twisting indicates that this disc will achieve rigid-body precession with a twist of . This is the first indication in our results that a disc which is predicted to be disrupted can nonetheless generate internal stresses to enforce rigid precession when the degree of twist becomes large.

Models with h=0.01 and R=10 showed different behaviour, largely as a result of the larger precession rate which is induced for a larger disc, and suggest that the models run with R=8 had outer radii slightly smaller than the natural tidal truncation radius. The disc with , which supports bending waves, was observed to undergo modest tidal truncation and to break into two distinct disc parts which precessed independently of one another, an outcome which is similar to one obtained using an SPH simulation by Larwood et al. (1996) but for a disc with somewhat different parameters. A narrow rim at the edge of the disc, running between radii r= 7-10, detached from the interior disc and precessed independently at a faster rate. The inner disc part precessed at a rate similar to that observed in the R=8 case. Interestingly we estimate that for this larger disc model, such that the breaking of the disc is indeed expected.

The model with propagates warps diffusively, and was observed to differentially precess through 180 degrees, but maintained a very smooth twist profile throughout the simulation. This disc is clearly disrupted through strong differential precession, and the internal stresses are insufficient for the disc to achieve rigid-body precession. Over longer times this disc should become completely twisted up by differential precession such that it loses all semblance of a disc-like structure.

The results for the R=10 discs are interesting, as they indicate two different modes by which a disc may break or be disrupted. In the bending wave regime, our results suggest that disc disruption may occur via a discontinuous change in the precession rate at a particular radius, such that the disc breaks into two distinct pieces which precess independently. Within each separate disc-piece warp communication allows for rigid-body precession. The detachment of the disc outer rim may in part be influenced by the strong tidal interaction there, since non linear features and shocks may have the effect of reducing the efficacy of warp communication. In the diffusion regime, it appears that viscosity is able to smooth out any discontinuities, and the disc is disrupted through global differential precession of the disc annuli, with a smooth twist profile being maintained across the disc radius.

A significant result obtained for the h=0.01 models is that these highly twisted discs align with the binary orbit plane much more quickly than the thicker discs, which align on the viscous timescale. Although we find that numerical diffusion can cause an inclined disc to align with the equatorial plane of the computational domain, this occurs on much longer timescales than are observed. We tentatively suggest that alignment of highly twisted discs occurs on the precession time rather than the viscous time, causing rapid alignment of very thin discs.

We have not considered different binary mass ratios, instead we restricted our study to the case of an equal mass binary. As Eqs. (38) and (15) indicate, the free and mean precession rates scale linearly with the companion mass for a disc of fixed outer radius. Hence we would expect the differential precession timescale to be increased for less massive binary companions, and consequently the differential twist observed in our models should be reduced in such systems.

Our work has a number of astrophysical implications. Most T Tauri stars are members of binary or multiple systems, and in a number of these it is believed that the disc and orbit planes are misaligned. One particular system for which there is a resolved image of a disc which is misaligned from the binary plane is HK Tau (Stapelfeldt et al. 1998). In this system, the disc radius is estimated to be 105 AU, and the projected binary separation is approximately 3 times larger than this, suggesting that the disc may be tidally truncated and undergoing strong interaction with the inclined companion star. The images do not show any signs that the disc is warped, which is consistent with our h=0.05 models, whose aspect ratio is probably appropriate for T Tauri discs. Evidence for misaligned young binaries is also provided by observations of precessing jets in star forming regions (Terquem et al. 1999). In one particular example, Bally & Devine (1994) suggest that the jet which appears to be launched by the young stellar object HH43* in the L1641 molecular cloud in Orion precesses with a period of about 10⁴ years. Applying models to a disc with outer edge of 50 AU surrounding a solar type star, the rigidly precessing disc models presented in this work precess with a period of between years, depending on the inclination angle considered. Thus the observed precession is consistent with that driven by a binary companion, with parameters close to those adopted in this work. The long alignment timescales that we find for h=0.05 discs suggests that a T Tauri disc will remain misaligned throughout its lifetime, such that jets launched from the disc inner regions will continue to precess for the duration of the time when jet launching is active.

The thin h=0.01 disc models that we have computed are most likely relevant to discs in X-ray binaries. The X-ray binary systems SS433 (Margon 1984) and Her X-1 (Boynton et al. 1980) are two examples of systems thought to contain precessing discs, with the evidence provided by the precessing jet of SS433 being particularly compelling. It is likely that the transfer of matter between the donor star and compact object in these systems will cause disc material to lie in the binary orbit plane initially, and a tilting mechanism is required to move the disc out of this plane, such as the Bardeen-Petterson effect (Bardeen & Petterson 1975), a misaligned dipole magnetic field associated with the central neutron star (Terquem & Papaloizou 2000), a disc wind (Schandl & Meyer 1994) or radiation pressure (Iping & Petterson 1990; Pringle 1996). If our results on the rapid realignment of highly twisted discs are correct, then this implies that a strong tilting mechanism which operates on short timescales will be required to maintain a misaligned disc which can be observed to precess.

Finally, we briefly consider the implications of our results for the formation of planets in misaligned binary systems. The early stages of planet formation are believed to involve the growth of planetesimals via mutual collisions. As the planetesimals grow in size, their interaction with the disc through gas drag forces will decrease. Planetesimals orbiting at different radii will have their own precession frequencies, likely to be different to that of the disc. Consequently, the growth of the planetesimals will eventually cause them to develop orbits which are increasingly inclined to the disc midplane. The effect of this is likely to be an increase in relative motion between planetesimals, which will affect the outcomes of mutual collisions, and the rate at which planetary growth proceeds. A detailed examination of these issues will be presented in a forthcoming paper.

References

All Tables

Table 1:   Table of runs.

All Figures

---

Copyright ESO 2010