Generalised warped disk equations

TL;DR

This paper introduces a unified set of one-dimensional equations that describe the evolution of warped accretion disks across both low and high viscosity regimes, addressing a gap in existing models.

Contribution

It provides a novel, simplified mathematical framework applicable to both wave-like and diffusive warp evolutions in accretion disks.

Findings

The equations work for arbitrary tilts with small warps.

They unify the description of warp evolution across viscosity regimes.

Potential for improved modeling of accretion disk dynamics.

Abstract

The manner in which warps in accretion disks evolve depends on the magnitude of the viscosity. For small viscosity $(\alpha < H/R)$ the warp evolves in a wave-like manner; for large viscosity $H/R<\alpha \ll 1$ it evolves diffusively. Here $\alpha$ is the viscosity parameter and $H/R$ the disk aspect ratio. Currently there is no simple set of equations which describes the evolution in both regimes. In this paper we describe a possible solution to this problem and introduce a set of one-dimensional equations that describe the evolution of a warped disk that are applicable in both high and low viscosity regimes for arbitrary tilts, but small warps.

Equations

$$\mathbf{L}=\Sigma R^2\Omega \mathbf{l},$$

$${\nu }_1={\alpha }_1{\left(\frac{H}{R}\right)}^2{R}^2\Omega$$

$${\nu }_2={\alpha }_2{\left(\frac{H}{R}\right)}^2{R}^2\Omega ,$$

$$\frac{\partial \mathbf{G}}{\partial t}+\omega \mathbf{l}\times \mathbf{G}+\alpha \Omega \mathbf{G}=\frac{\Sigma H^2{R}^3{\Omega }^3}{4}\frac{\partial \mathbf{l}}{\partial R}$$

$$\Sigma R^2\Omega \frac{\partial \mathbf{l}}{\partial t}=\frac{1}{R}\frac{\partial \mathbf{G}}{\partial R}+\mathbf{T},$$

$$\omega =\frac{{\Omega }^2-{\kappa }^2}{2\Omega }$$

$$\frac{\partial \mathbf{L}}{\partial t}=$$

$$+\frac{1}{R}\frac{\partial }{\partial R}\left[\left({\nu }_2{R}^2{\left|\frac{\partial \mathbf{l}}{\partial R}\right|}^2-\frac{3}{2}{\nu }_1\right)\mathbf{L}\right]$$

$$+\frac{1}{R}\frac{\partial }{\partial R}\left[\frac{1}{2}{\nu }_{2}R|\mathbf{L}|\frac{\partial \mathbf{l}}{\partial R}\right]+\mathbf{T}$$

$$\frac{\partial \Sigma }{\partial t}+\frac{1}{R}\frac{\partial }{\partial R}(R\Sigma v_R)=0,$$

$$\frac{\partial \mathbf{L}}{\partial t}+\frac{1}{R}\frac{\partial }{\partial R}(\Sigma v_R{R}^3\Omega \mathbf{l})=\frac{1}{R}\frac{\partial \mathbf{G}}{\partial R}+\mathbf{T},$$

$$v_R=\frac{\partial \mathbf{G}/\partial R\cdot \mathbf{l}}{R\Sigma d(R^2\Omega )/dR}.$$

$$\frac{\partial \Sigma }{\partial t}=-\frac{1}{R}\frac{\partial }{\partial R}\left[\frac{\partial \mathbf{G}/\partial R\cdot \mathbf{l}}{d(R^2\Omega )/dR}\right].$$

$$\frac{\partial \mathbf{L}}{\partial t}=$$

$$\frac{\partial \Sigma }{\partial t}=-\frac{2}{R}\frac{\partial }{\partial R}\left[\frac{(\partial \mathbf{G}/\partial R\cdot \mathbf{l})}{R\Omega }\right]$$

$$\frac{\partial \mathbf{L}}{\partial t}=-\frac{2}{R}\frac{\partial }{\partial R}\left[\left(\frac{(\partial \mathbf{G}/\partial R\cdot \mathbf{l})}{\Sigma R\Omega }\right)\mathbf{L}\right]+\frac{1}{R}\frac{\partial \mathbf{G}}{\partial R}+\mathbf{T}.$$

$$\frac{\partial \mathbf{G}}{\partial t}+\omega \mathbf{l}\times \mathbf{G}$$

$$\frac{\Sigma H^2{R}^3{\Omega }^3}{4}\frac{\partial \mathbf{l}}{\partial R}-\frac{3}{2}(\alpha +\beta ){\nu }_1\Sigma R^2{\Omega }^2\mathbf{l}.$$

$$\alpha \Omega \mathbf{G}+$$

$$\frac{\Sigma H^2{R}^3{\Omega }^3}{4}\frac{\partial \mathbf{l}}{\partial R}-\frac{3}{2}(\alpha +\beta ){\nu }_1\Sigma R^2{\Omega }^2\mathbf{l}.$$

$$\mathbf{G}=\frac{1}{2}{\nu }_2\Sigma R^3\Omega \frac{\partial \mathbf{l}}{\partial R}-\frac{3}{2}{\nu }_1\Sigma R^2\Omega \mathbf{l}.$$

$$\Sigma (R,0)={\Sigma }_0{\left(\frac{R}{R_{\mathrm{in}}}\right)}^{-1/2}\left[1-{\left(\frac{R_{\mathrm{in}}}{R}\right)}^{\frac{1}{2}}\right]\left[1-e^{R-R_{\mathrm{out}}}\right].$$

$$i(R,0)=10^{\circ }\left[\frac{1}{2}\tanh \left(\frac{R-R_{\mathrm{warp}}}{R_{\mathrm{width}}}\right)+\frac{1}{2}\right].$$

Full text

1. Introduction

Warped disks are expected to occur in a large number of astrophysical situations (e.g. Pringle, 1981, 1999; King et al., 2013). Warping may occur due to external torques from various sources. Binaries can provide such a torque on a circumstellar disk from an external binary component (e.g. Papaloizou & Terquem, 1995; Larwood et al., 1996; Lubow & Ogilvie, 2000; Martin et al., 2009, 2011) or a torque on a circumbinary disk form an internal binary (e.g. Facchini et al., 2013; Lodato & Facchini, 2013; Martin & Lubow, 2017). Around spinning black holes, general relativistic Lense–Thirring precession may cause warping (Bardeen & Petterson, 1975) in X–ray binaries (e.g. Scheuer & Feiler, 1996; Wijers & Pringle, 1999; Ogilvie, 2001; Martin et al., 2007) and around supermassive black holes (e.g. Herrnstein et al., 1996; Martin, 2008). Disks in AGN and in binary X-ray sources may be warped by the effects of radiation pressure (Pringle, 1996; Ogilvie & Dubus, 2001). Disks may also be warped by a misaligned magnetic field (e.g., Lai, 1999) or a planet (e.g. Nealon et al., 2018).

The evolution of the warp in a disk depends upon how the Shakura & Sunyaev (1973) viscosity $\alpha$ parameter compares with the disk aspect ratio $H/R$. If $\alpha>H/R$, then the warp propagates diffusively through the effects of viscosity. If $\alpha<H/R$, then pressure forces drive the evolution and the warp propagates as a bending wave that travels at half the sound speed, $c_{\rm s}/2$ (Papaloizou & Lin, 1995; Pringle, 1999; Lubow et al., 2002). In this case, the viscosity is too small to damp the wave locally. For a recent review of warped disks, see Nixon & King (2016).

One-dimensional (based on radius) models of disk tilt evolution offer advantages over multi-dimensional models. They permit tracking the evolution over long timescales with much less computational effort than multi-dimensional models. They are also easier to interpret physically. On the other hand, their applicability is limited by the simplifications made to reduce the dimensionality. In any case, such models can be compared with multi-dimensional models to obtain more physical insight.

A one-dimensional model should ideally conserve angular momentum and be valid for arbitrary tilts and warps (the derivative of the tilt with respect to the logarithm of the radius). In the viscous regime, Pringle (1992) developed a set of intuitively-based one-dimensional equations that satisfy these conditions. Ogilvie (1999) extended this analysis by directly working with the fluid equations and obtained one-dimensional equations that apply for even large warps. This analysis showed that the Pringle (1992) equations are valid for small warps and small $\alpha$, $H/r<\alpha\ll 1$, but arbitrary tilts, with the extension that the effective viscosity coefficients are constrained by the internal fluid dynamics.

In application to disks around young stars, the wave–like regime is of importance. One-dimensional linear disk evolution equations for this regime typically assume that the warp is small and ignore disk surface density evolution (e.g., Papaloizou & Lin, 1995; Lubow & Ogilvie, 2000). Ogilvie (2006) analyzed the nonlinear dynamics of free warps (imposed by initial conditions) in the absence of viscous density evolution. However, as found in Bate et al. (2000), significant density evolution can occur as the tilt evolves. We are interested in the case that the disk is in good radial communication so that the level of warping is small, as should apply to protostellar disks. Therefore a linear analysis is often valid. We are interested in the case that the disk tilt and surface density change over the course of its evolution.

The goal of this work is to find a formulation which describes the disk evolution correctly in both regimes, and manages to connect the two. In Section 2 we present two sets of warped disk equations, one valid in the viscous regime and one valid in the wave–like regime. In Section 3, we follow along the lines of Pringle (1992) to extend the linear tilt evolution equations to apply to arbitrary tilts and account for viscous density evolution. In Section 4 we numerically solve the equations for an initially warped disk around a single central object in the absence of any external torques. We draw our conclusions in Section 5.

2. Warped disk equations in the two regimes

Currently there are two sets of warped disk equations that describe mutually exclusive regimes, the wave–like regime with $\alpha<H/R$ and the diffusive (or viscous) regime with $\alpha>H/R$. We provide an overview of these two sets of equations in this Section. We describe the disk as consisting of a set circular rings with spherical radius $R$. We assume that the disk is in near Keplerian rotation. 111 Thus our results do not apply to strongly non–Keplerian flows such as those that occur in accretion disks close to the event horizon of a black hole. The rings rotate with Keplerian angular frequency $\Omega(R)=\sqrt{GM/R^{3}}$ about a central object of mass $M$ and have a surface density $\Sigma(R)$. The disk extends from inner radius $R_{\rm in}$ to outer radius $R_{\rm out}$. The angular momentum per unit area of each ring is

[TABLE]

where $\bm{l}$ is a unit vector. We consider a locally isothermal disk with aspect ratio $H/R$, where $H$ is the disk scale height.

Angular momentum transport in an accretion disk is driven by turbulent eddies with a maximum size $H$, and maximum speed the sound speed, $c_{\rm s}=H\Omega$. The azimuthal shear viscosity has the standard form

[TABLE]

(Shakura & Sunyaev, 1973), where $\alpha_{1}\simeq\alpha$, for dimensionless parameter $\alpha<1$. The vertical shear viscosity is

[TABLE]

(Papaloizou & Pringle, 1983), where $\alpha_{2}\simeq 1/(2\alpha)$ in the linear approximation.

2.1. Wave–like limit equations

In the wave–like limit, $\alpha<H/R$, we assume that the surface density does not evolve, $\partial\Sigma/\partial t=0$. The evolution of a warped disk is described by two equations

[TABLE]

and

[TABLE]

where $\bm{G}$ is the internal disk torque and $\bm{T}$ is the external torque on the disk (see equations 12 and 13 in Lubow & Ogilvie, 2000). The apsidal precession frequency in the plane of the disk is

[TABLE]

with epicyclic frequency $\kappa$. We solve equations (4) and (5) in Section 4.1 to compare to our solution to the generalised equations in the wave–like limit.

2.2. Viscous limit equations

In the viscous limit, $\alpha>H/R$, a warped disk is described by the evolution equation

[TABLE]

(Pringle, 1992). We solve equation (9) in Section 4.2 to compare to our solution to the generalised equations in the viscous regime.

3. Generalised warped disk equations

We now show how it is possible to combine the above equations into a single set that are valid in both the wave–like and the diffusive warp propagation regimes. We follow the methods of Papaloizou & Pringle (1983) and Pringle (1992). We note that this is not a first principles derivation of the evolution equations (cf., Ogilvie, 1999, 2006). Conservation of mass is expressed as

[TABLE]

where $v_{R}$ is the radial velocity. Conservation of angular momentum gives us

[TABLE]

where $\bm{G}$ is the internal disk torque and $\bm{T}$ is the external torque on the disk. In equation (11), we have included the second term on the LHS compared to equation (5), in the wave–like equations, in order to enforce conservation of angular momentum as the disk density evolves.

We take the dot product of equation (11) with $\bm{l}$ and subtract $R^{2}\Omega$ times equation (10) to obtain an equation for the radial velocity

[TABLE]

We substitute the radial velocity equation (12) into the conservation of mass equation (10) to obtain an equation for the surface density evolution

[TABLE]

Further, we substitute the radial velocity equation (12) into the conservation of angular momentum equation (11) to obtain an equation for the evolution of the angular momentum in the disk

[TABLE]

Since the disk is in near–Keplerian rotation we take $\Omega\propto R^{-3/2}$ and find

[TABLE]

and

[TABLE]

The key step to generalising the equations so that they are valid in both diffusive and wave–like regimes is to now amend equation (4) to read

[TABLE]

To effect this generalisation we have found it necessary to introduce two extra terms dependent on a new dimensionless parameter $\beta$. The fourth term on the left hand side has the effect of damping the component of disk torque $\bm{G}$ perpendicular to the local disk plane. The addition of the final term on the right hand side is to add an additional shear viscosity term. At this stage the magnitude of $\beta$ is arbitrary, except that we shall require that $\beta\gg\alpha$. We show the effects of different values for $\beta$ in Section 4.

3.1. The new generalised equation in the two limits

Equations (16), and (18) provide a one-dimensional description of both the the disk surface density and the disk tilt. We now show that this generalised equation has the previous equations (Section 2) in both limits.

In the wave-like limit we have $\alpha<H/R\ll 1$. The equations derived in this limit assumed that the surface density did not change with time, because in this limit the wave-like warp propagation happens on the shorter timescale than the viscous evolution of the surface density. Thus, in this limit, the final term on the right hand side of equation (18) is negligible. In addition the assumption that $\Sigma$ is independent of time implies that $v_{\rm R}=0$, unless there is an external source of mass. Thus (equation 12) we may take $\partial\bm{G}/\partial R\cdot\bm{l}=0$ and we may ignore the fourth term on the left hand side. Given this, equation (18) now reduces to equation (4), as required. We note, however, that the full solution to the new equations allows for the evolution of the surface density also in the wave-like regime. The degree to which the surface density evolves depends on the magnitude of the new parameter $\beta$. 2. 2.

In the viscous limit $(\alpha>H/R)$, $\bm{G}$ evolves on a viscous timescale and so $\partial\bm{G}/\partial t\ll\alpha\Omega\bm{G}$ and we set $\partial\bm{G}/\partial t=0$. Furthermore, we set $\omega=0$ provided that $\omega\ll\alpha\Omega$ and we are left with

[TABLE]

We take the dot product of this with $\bm{l}$ to find an expression for $\bm{G}\cdot\bm{l}$ and then substitute that into equation (20) to find

[TABLE]

Substituting this equation for $\bm{G}$ into equation (16), we recover the viscous disk evolution equation (9), which is valid for $H/r<\alpha\ll 1$ (Ogilvie, 1999).

4. Numerical solutions

We solve equations (16) and (18) as an initial value problem for $\bm{L}$, and $\bm{G}$ using finite differences. The method is first–order explicit in time. We use Cartesian coordinates and treat each component of the vectors separately. The units in the code are defined with $G=M=1$, where $M$ is the mass of the central object. The Keplerian orbital period at the inner disk radius $R_{\rm in}=1$ is $P_{\rm in}=2\pi$. We take the boundary conditions that $\bm{G}=\bm{0}$, $\Sigma=0$ and $\partial\bm{l}/\partial R=0$ at $R=R_{\rm in}$ and $R=R_{\rm out}$. The initial condition on $\bm{G}$ is always taken as $\bm{G}(R,0)=\bm{0}$.

We consider the evolution of an initially warped disk around a single central object. There is no external torque on the disk, so $\bm{T}=0$ and $\omega=0$. The disk extends from $R_{\rm in}=1$ up to $R_{\rm out}=20$. We take the initial surface density of the disk to be distributed as a simple power law with ends truncated at $R_{\rm in}$ and $R_{\rm out}$

[TABLE]

The constant $\Sigma_{0}$ is arbitrary as the equations are linear in $\Sigma$. Here we have scaled the total disk mass to be $0.001\,M$. The first factor in brackets on the RHS is a power law that represents a steady disk with $\nu_{1}\Sigma\propto\,\rm const$ if mass is added at the outer edge. The second and third factors enforce zero torque ($\Sigma=0$) inner and outer boundary conditions, respectively. Note that the surface density is not in steady state since we do not add material to the disk.

The initial tilt of the disk is described by

[TABLE]

Since the equations are linear in disk tilt, the normalisation of $i$ is arbitrary. The disk has an inclination of zero at the inner disk edge, an inclination of $10^{\circ}$ at the outer disk edge and a warp at radius $R_{\rm warp}=10$ with a width of radius $R_{\rm width}=2$. There is no twist in the disk. Since we do not have any torques to cause precession, the disk remains untwisted throughout its evolution. Thus we consider only the inclination of the disk and not the nodal precession angle.

4.1. Wave–like propagation; $\alpha<H/R$

We consider the evolution of an initially warped disk with parameters in the wave–like limit, $\alpha=0.01$ and $H/R=0.1$. Figure 1 shows the disk inclination and surface density evolution for several cases. In the top left panel we solve the wave–like warped disk equations (4) and (5) with a fixed density distribution. The warp in the disk propagates both inwards and outwards. The inwards propagating wave reflects off the inner boundary and then begins to propagate outwards.

In the other panels of Fig. 1 we solve the full disk equations (16) and (18) with different values for $\beta$. In the top right panel we show the behaviour that occurs if we do not introduce the parameter $\beta$. With $\beta=0$ the result is that there appears to be unphysical evolution of the disk surface density which occurs where the initial warp change was strongest, and which continues long after the initial warp has propagated away. The surface density anomaly shuld not keep growing at the position of the initial tilt change, even when the tilt at that point has evolved elsewhere. Furthermore, this behaviour is not seen in three dimensional hydrodynamical simulations (e.g. Nealon et al., 2015). This unphysical behaviour was the reason for introducing the new parameter $\beta$. The inclination evolution is very similar when we solve the wave–like equations (top left panel) or the full equations for $\beta\gtrsim 1$ (bottom panels). However, there is surface density evolution when we solve the full equations and angular momentum is conserved. The bottom two panels show that for $\beta\gtrsim 1$, the surface density evolution is independent of the value for $\beta$. There is slight difference between $\beta=1$ and $\beta=10$ but we find no difference for even higher $\beta$ compared to $\beta=10$.

4.2. Diffusive warp propagation; $\alpha>H/R$

As a check, we consider the evolution of a disk with parameters in the diffusive regime. We take $\alpha=0.1$ and $H/R=0.01$. Fig. 2 shows the disk inclination and surface density evolution solving the full equations (16) and (18) with $\beta=10$. We have also solved the diffusive equation (9) but find there is no difference between the two solutions and so we do not show this. In the diffusive regime there is no difference between solving the diffusive equations and the full equations that we have derived. The additional $\beta$ damping term has no effect in this limit.

4.3. Intermediate regime; $\alpha=H/R$

Neither the wave-like equations nor the diffusive equations are able to model the evolution of a disk with $\alpha\approx H/R$. However, the full equations we have developed, equations (16) and (18), can be used in this regime. Fig. 3 shows the solution to the full equations with $\beta=10$ for a disk with $\alpha=0.1$ and $H/R=0.1$. The inner parts of the disk appear more diffusive in nature and the outer parts look more wave-like in the inclination evolution.

5. Conclusions

We have introduced a new set of equations that describe the evolution of disk warp and of disk surface density in both low viscosity and high viscosity disks. We have shown that the two sets of equations agree with the equations for warp propagation previously derived in the two distinct regimes of low viscosity (wave-like warp propagation) and of high viscosity (diffusive warp propagation). In order to achieve this we have introduced a new dimensionless parameter $\beta$ which has the dominant effect of preventing unphysical evolution of surface density in the wave-like regime. We have not been able to determine the required magnitude of $\beta$ except to note that for $\beta\gg\alpha$ the unphysical evolution of surface density in the wave-like regime no longer occurs. In order to determine the value of $\beta$, and indeed to determine whether or not the new equations we present here provide an adequate description of warp evolution in general, it will be necessary to undertake a detailed analytic analysis (c.f. Ogilvie, 1999) and/or compare with detailed numerical simulations.

Acknowledgments

RGM, SHL and AF acknowledge support from NASA through grant NNX17AB96G. RN has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601). CJN is supported by the Science and Technology Facilities Council (grant number ST/M005917/1).