Long-Lived Eccentricities in Accretion Disks

TL;DR

This paper demonstrates that eccentric modes in accretion disks with realistic boundaries are trapped and long-lived, potentially explaining observed asymmetries in protoplanetary disks, through a simple WKB theoretical framework.

Contribution

It shows that in pressure-only accretion disks with realistic boundaries, eccentric modes are trapped and long-lived, providing a new understanding of disk asymmetries.

Findings

Eccentric modes are trapped in disks with rapid surface density drops.

The longest-lived mode is the zero-node aligned elliptical streamline mode.

Such modes decay on the viscous timescale, living for the disk's lifetime.

Abstract

Accretion disks can be eccentric: they support $m=1$ modes that are global and slowly precessing. But whether the modes remain trapped in the disk---and hence are long-lived---depends on conditions at the outer edge of the disk. Here we show that in disks with realistic boundaries, in which the surface density drops rapidly beyond a given radius, eccentric modes are trapped and hence long-lived. We focus on pressure-only disks around a central mass, and show how this result can be understood with the help of a simple second-order WKB theory. We show that the longest lived mode is the zero-node mode in which all of the disk's elliptical streamlines are aligned, and that this mode decays coherently on the viscous timescale of the disk. Hence such a mode, once excited, will live for the lifetime of the disk. It may be responsible for asymmetries seen in recent images of protoplanetary disks.

Equations

$$\omega E=\frac{1}{2\Omega r^3\Sigma }\left\{\frac{d}{dr}\left(\gamma r^{3}P\frac{dE}{dr}\right)+r^2\frac{dP}{dr}E\right\},$$

$$y=(r^{3}P{)}^{1/2}E,$$

$$\frac{d^{2}y}{d{r}^2}+\frac{2\Omega }{c^2}\left({\omega }_p(r)-\omega \right)y$$

$$c=\sqrt{\gamma P/\Sigma }$$

$${\omega }_p(r)=-\frac{c^2}{2\Omega }\left[(r^{3}P{)}^{-1/2}\frac{d^2}{d{r}^2}(r^{3}P{)}^{1/2}-\frac{1}{\gamma rP}\frac{dP}{dr}\right].$$

$$T(r)$$

$$\Sigma (r)$$

$$\omega ={\omega }_p-\frac{c^2}{2\Omega }{k}^2,$$

$$\text{varointclockwise}kdr=(2n+1)\pi =\pi ,3\pi ,5\pi ,\cdots ,$$

$${\omega }_0\equiv \frac{c_0^2}{2{\Omega }_0{r}_0^2},$$

$$dE/dr=0.$$

$$\omega =-\left(n+\frac{1}{2}\right)+{\omega }_{p,\mathrm{peak}},$$

$$(2n+1)\pi$$

$$\sim$$

$$|\omega |\sim n^{1/\delta },$$

$$E\sim (r^{3}P{)}^{-1/2}\exp \left[-S(r)\right],$$

$$E\sim r^{\frac{p+q-3}{2}-\frac{B}{2-p}}\exp \left[-\frac{\omega r^{p+q-3/2}}{(2-p)(p+q-3/2)}\right].$$

$$\displaystyle\Gamma=\frac{\int\alpha r^{3}P\big{|}dE/dr\big{|}^{2}dr}{\int 2\Omega r^{3}\Sigma|E|^{2}dr}\sim\frac{\int\alpha\big{|}dy/dr\big{|}^{2}dr}{\int(2\Omega/c^{2})|y|^{2}dr},$$

$${\omega }_p=-r^{-q-1/2}(A+B{r}^{\xi }+C{r}^{2\xi }).$$

$$A$$

$$B$$

$$C$$

Full text

Protoplanetary disks (1300), Circumstellar disks (235), Eccentricity (441), Astrophysical fluid dynamics (101)

1 Introduction

Accretion disks are usually assumed to be circular. But for disks in nearly Keplerian potentials, orbits are in general eccentric, and hence the disk as a whole might be eccentric. Eccentric disks are interesting for a variety of reasons; e.g., their distorted shape could be observed, and planets born within such disks would be eccentric, perhaps explaining observed planetary eccentricities.

Eccentric orbits of fluid around a star precess differentially due to the effects of pressure, self-gravity, and non-Keplerian potential components. However, particular eccentricity profiles may be found such that the disk as a whole precesses rigidly while maintaining its eccentricity. To find such profiles, one may linearize the equation of motion for the (complex) eccentricity, in which case the solution is a sequence of normal modes. The “fundamental” zero-node mode, i.e., the one in which the eccentricity is never zero, is typically the one of primary interest (e.g., it lives longest in the presence of viscosity).

The aforementioned normal mode calculation has been studied by many authors (Lubow, 1991; Tremaine, 2001; Papaloizou, 2002; Goodchild & Ogilvie, 2006; Ogilvie, 2008; Saini, Gulati, & Sridhar, 2009; Lubow, 2010; Teyssandier & Ogilvie, 2016). In Lee, Dempsey, & Lithwick (2019) (hereafter Paper I), we solved the normal mode problem in disks subject to both pressure and self-gravity forces, and also explained the numerical results in terms of a simple WKB theory. However, an important problem with virtually all previous studies is the boundary condition. It is usually assumed that the disk has sharp edges, beyond which the surface density drops instantaneously to zero. But real disks likely have a more gradual drop (Lynden-Bell & Pringle, 1974). One might worry that normal modes would not be trapped in such a disk, leading to no solutions in which the disk can remain eccentric.

In this paper, we shall show that this worry is unfounded. Virtually any disk with a realistic surface density profile will lead to trapped eccentric modes. Furthermore, we calculate with numerics and theory the eccentricity profile and precession rate that occur in disks of arbitrary surface density profiles. Throughout this paper, we focus on the pressure-only case, i.e., we ignore self-gravity. We do this for simplicity, but also because even if self-gravity is important in the bulk of the disk, in the outer part where the surface density is small pressure will play the larger role (Paper I). We organize this paper as follows. In Section 2, we present the equations of motion and their WKB formulation. In Section 3 we present numerical eigen-solutions, and show how these results may be understood with WKB theory. In Section 4 we address additional effects, before concluding in Section 5.

2 Equation of Motion

We consider a two-dimensional gas disk orbiting a central object of mass $M_{\star}$. The disk is cold and thin, such that its aspect ratio is much less than one. Perturbed variables are assumed to depend on time and azimuth in proportion to $e^{-i\omega t+m\varphi}$. We further set $m=1$, as is the case for eccentric perturbations, in which case $\omega$ may be identified as the precession frequency. Using the equations of continuity, momentum, and entropy, one can derive the following eccentricity equation by expressing the fluid variables in terms of eccentricity (Goodchild & Ogilvie 2006; Ogilvie 2008; Teyssandier & Ogilvie 2016; see Appendix A of Paper I):

[TABLE]

where $E$ is the complex eccentricity, $\Sigma$, $P$, and $\Omega=\sqrt{GM_{*}/r^{3}}$ are the surface density, two-dimensional pressure, and rotational frequency of the gas, respectively, and $\gamma$ is the adiabatic index. The complex eccentricity $E=|E|e^{-i\varpi}$ is used because both the amplitude $|E|$ and the periapse angle $\varpi$ are radial functions. We assume the perturbation is adiabatic111Corresponding equations for locally-isothermal perturbations (i.e., short cooling time) and vertically-integrated 3D disks can be found in Teyssandier & Ogilvie (2016)., while the equilibrium disk may have non-constant background entropy. In deriving Equation (1) we assume the mode frequency $\omega$ is much smaller than the orbital frequency $\Omega$ of the disk, i.e., the mode precesses very slowly.

2.1 Trapped Modes

Equation (1) may be cast into a more transparent form by transforming variables from $E$ to

[TABLE]

which leads to an equation with no single-derivative ($dy/dr$) term:

[TABLE]

where

[TABLE]

is the sound speed, and the “effective potential” $\omega_{p}$ is given by

[TABLE]

Equation (3) is particularly simple to analyze because of its similarity to Schrödinger’s equation, with $\omega_{p}$ playing the role of the potential and $\omega$ the role of the energy222 Unlike Schrödinger’s equation, Equation (3) has a spatially-variable coefficient ($2\Omega/c^{2}$) multiplying its energy term. That may be fixed by transforming co-ordinates (e.g., Ogilvie, 2008; Saini et al., 2009), but we do not do so.. It is apparent that $y$ is wave-like where $\omega_{p}>\omega$ and evanescent where $\omega_{p}<\omega$, with turning points at $\omega_{p}=\omega$.

For realistic disk profiles, the “potential” $\omega_{p}$ almost always has an inverted-U shape. As a result, eccentric modes are trapped by the peak of the potential, and their character will not depend on what happens far from the peak—a central result of this paper. For example, we plot $\omega_{p}$ in the middle panel of Figure 1 for our “fiducial case”: a disk with the temperature and surface density profiles $T\propto r^{-1/2}$ and $\Sigma=r^{-1}e^{-r}$ (Figure 1, top panel), respectively. More generally, motivated by the self-similar solutions of Lynden-Bell & Pringle (1974), we consider profiles of the form

[TABLE]

i.e., with a cutoff on the surface density at $r\gtrsim r_{0}$, and two free parameters, $p$ and $q$. For simplicity, we set $\xi$ in the above equations to $\xi=2-p$, as is the case for a self-similarly evolving disk that has power-law viscosity333If the disk also has an $\alpha$-viscosity with constant $\alpha$, then $q=3/2-p$. . The resulting $\omega_{p}$ is a sum of three power-laws (given explicitly in Equation (A1)). In Figure 2 the white region shows where, in the $p$-$q$ plane, $\omega_{p}$ has an inverted-U shape, and so can trap modes. In the grey region $\omega_{p}$ rises continually outwards, and so modes cannot be trapped. For the remainder of this paper, we consider only the white region because it encompasses most typically assumed values for $p$ and $q$.

2.2 Second-order WKB Theory

An advantage of analyzing Equation (3) rather than (1) is that one may derive a second-order-accurate WKB dispersion relation by replacing $d/dr\rightarrow ik$ (Gough, 2007):

[TABLE]

where $k$ is the radial wavenumber. We consider Equation (7) as a second-order dispersion relation because $\omega_{p}$ is two orders of $kr$ smaller than the leading term (Paper I).444Equation (1) admits a conserved, i.e., spatially constant, quantity $F=-i{m\pi\gamma\over 2}(y^{*}{dy\over dr}-y{dy^{*}\over dr})$, that is equal to the angular momentum flux. In WKB, $F=m\pi\gamma r^{3}Pk|E|^{2}$ (Goldreich & Tremaine, 1979).

Following Paper I, the dispersion relation can be analyzed using a frequency diagram (Figure 1, middle panel) and dispersion relation map (DRM; Figure 1, bottom panel), which plots contours of constant $\omega$ in the $r$-$kr$ plane. The frequency diagram shows the trapped modes. The mode with highest frequency is the zero-node fundamental mode, and those with more negative frequency have increasing number of nodes. Note that in this example (and as is typically true in the absence of self-gravity forces) frequencies are negative, meaning modes are retrograde. From the DRM, we see how a trapped wave refracts as it propagates, transitioning from outwardly-propagating (in the lower half of the figure) to inwardly propagating (top half) at its turning point ($k=0$). Its group velocity is $v_{g}=(\partial\omega/\partial k)_{r}=-kc^{2}/\Omega$.

In order for a trapped wave to represent a standing mode, its complex phase must change by an integral multiple of $2\pi$ over the course of a loop in the DRM. Accounting for the phase change at turning points, that implies that standing modes must satisfy

[TABLE]

where the path-integral is taken in the clockwise direction in phase-space (according to the group velocity), and $n=0,1,2,\cdots$ labels the number of nodes (Mark 1977; Shu et al. 1990; Paper I). Equation (8) allows us to determine the mode frequency from theory, after inserting $k(r;\omega)$ from the dispersion relation.

3 Solutions

3.1 Numerical Solutions

We solve Equation (1) numerically; in subsequent subsections we shall explain the results with second-order WKB theory. For the background state, we adopt Equations (5)–(6), with $p=1$, $\xi=1$, and $q=1/2$, and explore other values for these parameters below. We set the 2D adiabatic index to $\gamma=3/2$. Using $P=\Sigma c^{2}/\gamma$ in Equation (1) implies that its right-hand side is proportional to

[TABLE]

where $c_{0}$ and $\Omega_{0}$ are evaluated at the density cutoff $r_{0}$. The quantity $\omega_{0}$ is the characteristic frequency that determines the precession frequency of the modes, aside from dimensionless constants. We can write $\omega_{0}\sim h_{0}^{2}\Omega_{0}$ where $h_{0}$ is the aspect ratio of the disk at the outer cut-off, i.e., it is slower than the orbital frequency at the disk’s outer cutoff by the square of the aspect ratio there. Henceforth we shall measure mode frequencies in units of $\omega_{0}$ (i.e., set $\omega_{0}\rightarrow 1$). We choose our length unit to be $r_{0}$, i.e., set $r_{0}\rightarrow 1$.

For boundary conditions, the inner and outer disk edges are assumed to be free surfaces (i.e., zero Lagrangian pressure), such that the following boundary condition holds (Ogilvie 2008, Teyssandier & Ogilvie 2016, Paper I):

[TABLE]

The inner and outer radii of the disk are $r_{\rm min}=10^{-4}$ and $r_{\rm max}=50$, respectively.

We solve the boundary-eigenvalue problem of Equations (1) and (10) with a finite difference method (see Paper I for details) and a Chebyshev spectral method (Trefethen, 2000). In both cases, we construct a square matrix and solve for the eigenvalues and eigenfunctions using the LAPACK library. The two solution methods give essentially the same results, and we have checked that the results have converged with respect to the number of grid points.

3.2 Eigenfunctions

Figure 3 shows the numerical eigenfunctions for the first few modes. From the right-hand panel, which shows the scaled eccentricity $y$ (Equation (2)), we see that the oscillatory and evanescent regions predicted from the frequency diagram (Figure 1) match up with those of the solutions. The left-hand panel shows the unscaled eccentricity $E$, which rises outwards, as explained in further detail in Section 3.4.

3.3 Eigenvalues

In Figure 4, the orange dots show the numerically calculated eigenfrequencies up to $n=15$. Also shown, as blue dots, is the result from reducing the outer boundary from $r_{\rm max}=50$ to 10. For small $n$ the two sets agree, while for larger $n$ they diverge. The reason for the divergence is that, in the lower $r_{\rm max}$ simulation, the turning points no longer fit in the simulation domain when $n\gtrsim 3$ (Figure 3, right panel).

We may calculate the eigenfrequencies theoretically by inserting $k(r;\omega)$ from the dispersion relation (Equation 7) into the quantum condition (Equation 8). For the fiducial case $(p,\xi,q)=(1,1,1/2)$ (red circle in Figure 2), the integral can be performed analytically, yielding a harmonic-oscillator-like form:

[TABLE]

where $\omega_{p,\rm peak}=-(B+2\sqrt{AC})\simeq-0.81$ is the peak value of $\omega_{p}$. Figure 4 shows that the WKB formula above provides an excellent match to the numerical eigenvalues. In fact, it might appear that the agreement is too good at small $n$, given that we are making the WKB approximation, which one might expect to fail for the low $n$ modes. The reason is that, near the peak of $\omega_{p}$, the exact equation (Equation 3) is similar to the equation for a harmonic oscillator, for which the WKB solution is exact.

We now consider the more general case where $p$ and $q$ are arbitrary (but $\xi=2-p$), proceeding approximately. The quantum condition reads

[TABLE]

where $r_{\pm}$ are the turning points, and in the latter expression we dropped $\omega_{p}$ because it becomes small (in magnitude) relative to $\omega$ far from the turning points. To estimate $r_{+}$ we set $\omega=\omega_{p}$, where for $\omega_{p}$ we use its dominant piece, i.e., the last term in Equation (A1). We find

[TABLE]

where $\delta=(2q+1)/(14-8p-4q)+1/2$ and we omit order-unity coefficients.

3.4 Behavior of Eigenfunction at large radii

Figure 3 shows that at $r\gtrsim 1$ the eccentricity continues to rise outwards. Since the behavior near $r\sim 1$ is potentially observable, we examine it here in more detail. Figure 5 shows the fundamental modes for a set of background profiles with different $p$’s. Sometimes the eccentricity is peaked within the disk, while sometimes it rises continuously outwards. For more general values of $p$ and $q$, we have determined numerically where the eccentricity of the fundamental mode peaks. Our results are shown as blue dashed lines in Figure 2, labeled by $r_{\rm peak}$ (the radius where they peak). Curiously, beyond the constant $\alpha$ line the eccentricity peaks at infinite $r$; otherwise, it peaks further in.

We may understand this behavior from WKB theory. The solution beyond the outer turning point is a decaying exponential function for $y$ (e.g., Bender & Orszag, 1999), which leads to

[TABLE]

where $S=\int^{r}_{r_{+}}\sqrt{-k^{2}(s)}ds$ and we assume $r_{+}\ll r$. Inserting the expression for $k$ from the dispersion relation, and Taylor expanding in small $1/r$ yields $\sqrt{-k^{2}}\approx C^{1/2}r^{\xi-1}\left[1+{1\over 2C}(Br^{-\xi}+\omega r^{q-2\xi+1/2})\right]$. Note that the last term is subdominant because $q-2\xi+1/2<0$ in the white “trapped mode” zone in Figure 2. Therefore, we get

[TABLE]

We see that for $p+q-3/2>0$ (i.e., beyond the constant $\alpha$ line), the eccentricity rises inexorably outwards, while in the opposite limit it is peaked in the disk, confirming the numerical result.

4 Discussion

4.1 Viscous Damping

To study eccentricity damping, we consider the effects of bulk viscosity, following Goodchild & Ogilvie (2006). We ignore shear viscosity because it may cause overstability (Ogilvie, 2001). The linearized equation of motion with bulk viscosity is obtained by replacing the adiabatic index $\gamma$ by $\gamma+i\alpha$ in Equation (1), where $\alpha$ is the Shakura-Sunyaev parameter (but for bulk viscosity rather than the more common shear viscosity). We solve this modified equation numerically in the same way as before, and obtain the damping rate as the imaginary part of the eigenfrequency $\Gamma=-{\rm Im}\,\omega$.

We find that the numerical damping rate is of order of $\alpha|\omega|$. This can be understood by considering an integral relation from Equation (1) (Goodchild & Ogilvie, 2006; Lubow, 2010; Teyssandier & Ogilvie, 2016):

[TABLE]

where the absolute sign denotes the amplitudes. The last approximation is justified because the numerator is dominated by the oscillatory region within the wave cavity (Figure 3). Under the WKB approximation, the integrals can be estimated as $\Gamma\sim\alpha k^{2}c^{2}/2\Omega\sim\alpha|\omega|$, consistent with the numerical result.

This damping rate $(\Gamma\sim\alpha|\omega|)$ is slow. In particular, since $\omega\sim\omega_{0}$, the damping time $\Gamma^{-1}$ is comparable to the viscous time of the disk as a whole. Therefore, once such a mode is excited, it will live for the lifetime of the disk.

4.2 Three-dimensional Effects

We have assumed throughout this paper that the disk is two-dimensional. But Ogilvie (2008) pointed out that in 3D disks an extra precession term needs to be included because a fluid parcel cannot maintain vertical hydrostatic equilibrium while in an eccentric orbit. He derived an equation analogous to Equation (1) valid for 3D disks. We have analyzed that equation in the same way as we analyzed Equation (1), but we do not present details here—primarily because in order to derive a reduced 1D equation (such as Equation (1)), one must make the questionable assumption that the eccentricity is independent of height. Nonetheless, the result is very similar to the one found in this paper, i.e., there is an $\omega_{p}$ function that traps the modes. But there is one important difference, which is that in the inner disk the sign of the dominant term of $\omega_{p}$ switches to positive, which can potentially remove the inner turning point (the outer turning point is unaffected), and hence the modes must rely on the reflection at the disk inner edge in order to remain trapped (e.g., Miranda & Rafikov, 2018).

5 Conclusion

We demonstrate that typical accretion disks with a realistic outer density drop support trapped eccentric modes. We use the second-order WKB theory developed in Paper I to explain the different features, such as the wave cavity and the eccentricity in the outer part of the disk. 2. 2.

Each normal mode solution corresponds to a rigidly precessing eccentric pattern. This is a balance between the pressure precession effects caused by the axisymmetric and non-axisymmetric components of the disk. We find that, instead of the “test-particle precession rate” (i.e., the difference between orbital and epicyclic frequencies $\dot{\varpi}=\Omega-\kappa$) in celestial mechanics, the $\omega_{p}$ function gives the correct behavior of eccentricity of a gas disk. Previous results based on the leading-order WKB dispersion relation (e.g., Papaloizou, 2002; Goldreich & Sari, 2003) may require extra examination. 3. 3.

We find that trapped eccentric modes are standing waves with a discrete spectrum (Figure 1). The fundamental mode has zero radial nodes and the least negative frequency. It has the slowest damping rate when viscous damping is considered (Section 4.1). 4. 4.

The eccentricity of a mode in the outer disk is evanescent and does not carry angular momentum. It can also be explained by the WKB theory (Section 3.4). The eccentricity of a disk with $p+q<3/2$ peaks inside the disk rather than at infinity (Figure 2). 5. 5.

The trapped modes are not affected by the boundary conditions as long as the artificial computational boundary (i.e., $r_{\rm max}$) is far from the turning point (Figure 4).

Although we have shown that eccentric modes can live for a long time, we have not addressed the question of how such modes are excited. Some possibilities are gravitational excitation by a planet or star, or an internal instability.

Finally, we note that a long-lived eccentric disk can possibly be detected directly. Some lopsided disks have been detected (e.g., Dong et al., 2018; Andrews et al., 2018), and a variety of mechanisms proposed to explain them (e.g., Hsieh & Gu, 2012; Ataiee et al., 2013; Lyra & Lin, 2013; Mittal & Chiang, 2015; Zhu & Baruteau, 2016). But whether the lopsidedness might be due to the disk being eccentric is an intriguing possibility.

WKL thanks Kenny L.S. Yip for checking Equation (11) using contour integration. Y.L. acknowledges NASA grant NNX14AD21G and NSF grant AST-1352369.

Appendix A $\omega_{p}$ in model disk

For a disk model given by Equations (5)–(6) with general parameters $(p,q,\xi)$, the $\omega_{p}$ function in Equation (4) is given by

[TABLE]

The coefficients are given by

[TABLE]

when units are set by $\omega_{0}=1$ and $r_{0}=1$. For the fiducial case ($p=2-\xi=1$, $q=1/2$, and $\gamma=3/2$), $A=13/16$, $B=-1/12$, and $C=1/4$.