IRRADIATION INSTABILITY AT THE INNER EDGES OF ACCRETION DISKS

Jeffrey Fung; Pawel Artymowicz

doi:10.1088/0004-637X/790/1/78

1. INTRODUCTION

Accretion disks are susceptible to a wide range of instabilities, including the magnetorotational instability (Balbus & Hawley 1998), gravitational instability (Lin & Pringle 1987; Gammie 2001), Papaloizou–Pringle instability (Papaloizou & Pringle 1984, 1985, 1987; Goldreich et al. 1986), and Rossby wave instability (RWI; Lovelace et al. 1999). The list goes on as non-ideal MHD and vertical shearing (Urpin & Brandenburg 1998) are considered. These instabilities drive the evolution of disks by generating turbulence and creating complex, sometimes extreme structures, such as the formation of planets in protoplanetary disks.

Radiation pressure is a force generally present in all types of accretion disks. Its effect on accretion disks has been studied in many different aspects, including driving disk winds in active galactic nuclei (AGNs; e.g., Higginbottom et al. 2014), shaping particle size distributions in debris disks (Thebault et al. 2014), and influencing the motions of the inner rims of transitional disks (Chiang & Murray-Clay 2007; Dominik & Dullemond 2011). We demonstrate in this paper that radiation pressure can also cause a disk instability of its own kind. In the following, we give a brief introduction to this instability before launching into the formal theoretical work.

The strength of radiation pressure compared to gravity is measured by the number β:

$\begin{equation} \beta = \frac{\kappa _{\rm opa} L}{4 \pi cGM}, \end{equation} \tag{ 1 }$

where L is the central object's luminosity, M is its mass, and κ_opa is the opacity of the disk material; c and G are the speed of light and gravitational constant, respectively. The key to this instability is shadowing. As the front part of the disk gets pushed by radiation pressure, it also casts a shadow that reduces the amount of radiation pressure on the material behind the disk. In a one-dimensional (1D) radial picture, since radiation pressure always diminishes outward, the inner part of a disk always feels a stronger push than the outer part, and the net effect is therefore radial compression. In other words, any two concentric disk annuli would feel an attraction between them due to the combined effects of radiation pressure and shadowing.

This 1D scenario does not easily extend to a two-dimensional (2D) disk, however, because radiation pressure from a central source does not exert any azimuthal force. By the conservation of angular momentum, when a disk element is perturbed radially, it will oscillate at some epicyclic frequency. Figure 1 illustrates what effect this oscillating element has on the disk. One can see that disk material near the orbit of the perturbed element will experience a variation in shadowing along the azimuth. This variation generates a forcing that induces the unperturbed material to follow the motion of the perturbed element. The result is a global collective motion that is capable of growing on its own. We term this phenomenon the "irradiation instability" (IRI), since it relies on irradiation by the central object.

**Figure 1.** Simple illustration describing IRI. The blue curve denotes the orbit of a perturbed disk element oscillating around its guiding center, denoted by the dashed black line at r₀. The shaded area is where the disk sees the shadow cast by the perturbed element. The red arrows show the directions of radial forcing on the background disk relative to the average amount of radiation pressure received along r₀. These arrows are inward when they are in the shadow of the element, and outward when they are not. One can see that the background disk near r₀ is forced in the direction of amplifying the initial perturbation.
Download figure:
Standard image High-resolution image

Because a larger β allows for a more rapid radial motion, its value is crucial for the survival of this collective motion against disk shear. In most systems, dust grains provide the largest contribution to β. In circumstellar disks, micron-size grains can have β > 1 for F-type stars, and up to β ∼ 10¹ for A-type stars (e.g., Equation (10) of Kirchschlager & Wolf 2013). Given that the gas-to-dust ratio is typically ∼10², β of a perfectly coupled gas+dust mixture may be of the order of a few percent. Additionally, dust settling can enhance β in the midplane by reducing the local gas-to-dust ratio, while the radial migration of dust results in size segregation (Thebault et al. 2014), which can also enhance β at local radii. In other systems where radiation pressure can drive significant mass loss, such as AGN accretion disks, one would even expect β to exceed unity.

This paper aims to provide a basic understanding of IRI, of both the conditions that trigger it, and its consequences. In Section 2, we present a theoretical foundation for IRI and derive its instability criterion. Section 3 contains our disk model. Section 4 describes our semi-analytic and numerical methods. Section 5 reports the modal growth rate as a function of β and the sound speed c_s of the disk, and gives a discussion on the nonlinear evolution of IRI. Section 6 concludes with an outlook for future work.

2. THE LINEAR THEORY

We follow the method of Goldreich & Tremaine (1979), using similar notation, to derive the linear response of a 2D disk stirred by radiation pressure. We start with the continuity equation and the conservation of momentum:

$\begin{equation} \frac{\partial \Sigma }{\partial t} + {\bf \nabla }\cdot (\Sigma {\bf v}) = 0, \end{equation} \tag{ 2 }$

$\begin{equation} \frac{\partial {\bf v}}{\partial t} + ({\bf v}\cdot {\bf \nabla }){\bf v} = -{\bf \nabla }\eta -\frac{GM(1-\beta e^{-\tau })}{r^2}\hat{\bf r}, \end{equation} \tag{ 3 }$

where Σ is the surface density of the disk; v is the 2D velocity field; η is the specific enthalpy such that ∇η = ∇P/Σ, where P is the vertically averaged gas pressure; and τ is the optical depth of the disk. We denote the Keplerian orbital frequency as Ω_k, and the sound speed c_s is defined by the ideal gas law $P=c_{ s}^2 \Sigma$ . τ depends on the density distribution by the following equation:

$\begin{equation} \tau = \int ^r_0 \kappa _{\rm opa}\rho dr^{\prime }, \end{equation} \tag{ 4 }$

where ρ is the density of the disk. Near the midplane, ρ∝Σ/h, where h = c_s/Ω_k is the scale height of the disk. Note that with Equation (3) we have neglected the scattering of light into the azimuthal direction.

Σ, η, and $\bf v$ can be separated into a background quantity (without any subscript) and a perturbed quantity (denoted by the subscript "m"). We assume the background disk to be axisymmetric and in hydrostatic equilibrium so that v =(0, rΩ), where

$\begin{equation} \Omega = \sqrt{\Omega _{ k}^2 (1-\beta e^{-\tau })+\frac{1}{r}\frac{d\eta }{dr}}, \end{equation} \tag{ 5 }$

and the components of the perturbed velocity are denoted as v_m ≡( $u,\upsilon$ ). To simplify the notation, we also define the background and perturbed radiation force as F and F_m:

$\begin{equation} F = r \Omega _{ k}^2 \beta e^{-\tau }, \end{equation} \tag{ 6 }$

$\begin{equation} {F}_{m} = -F\int ^r_0 \frac{{\eta }_{m}}{c_{ s}^2}\frac{d\tau }{dr^{\prime }}dr^{\prime }. \end{equation} \tag{ 7 }$

For a small perturbation, it follows from Equations (2) and (3) that the perturbed quantities are governed by the following linearized equations:

$\begin{equation} \frac{\partial {\Sigma }_{m}}{\partial t} + {\bf \nabla }(\Sigma {\bf v}_{m}) + {\bf \nabla }({\Sigma }_{m}{\bf v}) = 0, \end{equation} \tag{ 8 }$

$\begin{equation} \frac{\partial {\bf v}_{m}}{\partial t} + ({\bf v}\cdot {\bf \nabla }){\bf v}_{m} + ({\bf v}_{m}\cdot {\bf \nabla }){\bf v} = -{\bf \nabla }{\eta }_{m} + {F}_{m} \hat{\bf r}. \end{equation} \tag{ 9 }$

Without a loss of generality, we can assume a form of the solution for the perturbed quantities Σ_m, η_m, u, and $\upsilon$ :

$\begin{equation} {X}_{m}(r,\theta,t)=X(r)e^{i\left(m\theta -\omega t\right)}, \end{equation} \tag{ 10 }$

for some complex function X(r) and complex number ω, while m is the azimuthal mode number. Substituting this form into Equation (9), we find

$\begin{equation} u = -\frac{i}{D}\left[\frac{2m\Omega }{r}{\eta }_{m} + {\Omega }_{m}\left(\frac{\partial {\eta }_{m}}{\partial r}-{F}_{m}\right)\right], \end{equation} \tag{ 11 }$

$\begin{equation} \upsilon = \frac{1}{D}\left[\frac{m{\Omega }_{m}}{r}{\eta }_{m} + 2\left(\Omega +\frac{r}{2}\frac{d\Omega }{dr}\right)\left(\frac{\partial {\eta }_{m}}{\partial r}-{F}_{m}\right)\right]. \end{equation} \tag{ 12 }$

The pattern rotation frequency Ω_m and the coefficient D are defined as

$\begin{equation} {\Omega }_{m}\equiv m\Omega -\omega, \end{equation} \tag{ 13 }$

$\begin{equation} D\equiv \kappa ^2-{\Omega }_{m}^2, \end{equation} \tag{ 14 }$

$\begin{equation} \kappa ^2= \frac{1}{r^3}\frac{d[r^4 \Omega ^2]}{dr}, \end{equation} \tag{ 15 }$

where κ is the epicyclic frequency of the unperturbed orbit. To solve for Σ_m, or equivalently η_m, we substitute Equations (11) and (12) into Equation (8), giving

$\begin{equation} \frac{\partial ^2{\eta }_{m}}{\partial r^2}+a(r)\frac{\partial {\eta }_{m}}{\partial r}+b(r){\eta }_{m}+c(r)\int ^r_0 \frac{{\eta }_{m}}{c_{ s}^2}\frac{d\tau }{dr^{\prime }} dr^{\prime } = 0, \end{equation} \tag{ 16 }$

where

$\begin{eqnarray*} \nonumber a &\equiv \frac{\partial }{\partial r}\ln {\left(\frac{r\Sigma }{D}\right)},\\ \nonumber b &\equiv \frac{2m\Omega }{r {\Omega }_{m}}\frac{\partial }{\partial r}\ln {\left(\frac{\Sigma \Omega _0}{D}\right)}-\frac{m^2}{r^2}+\frac{1}{c_{ s}^2}\left(F\frac{d\tau }{dr}-D\right),\\ \nonumber c &\equiv F \left(\frac{\partial }{\partial r}\ln {\left(\frac{r\Sigma F}{D}\right)}-\frac{2m\Omega }{r{\Omega }_{m}}\right). \end{eqnarray*}$

We arrive at a second-order integro-differential equation for η_m.

2.1. Instability Criterion

A local criterion for axisymmetric instability can be derived from Equation (16). We apply the WKB approximation and write ${\eta }_{m}\sim e^{i\int ^r_0 k_{ r} dr^{\prime }}$ , where k_r ≫ 1/r is the radial wave number. We then separate the real and imaginary part of the equation. Finally, setting m = 0, the dispersion relation can be written as

$\begin{equation} \omega ^2 = \kappa ^2 + k_{ r}^2 c_{ s}^2 - \Omega _{ k}^2 \beta e^{-\tau } \left(\frac{d\tau }{d\ln {r}} + {\tilde{\tau }}_{m}\frac{d\ln {[r{\mathscr R}]}}{d\ln {r}} \right), \end{equation} \tag{ 17 }$

where

$\begin{equation} \mathscr{R} \equiv \frac{\Sigma \Omega _k \beta e^{-\tau }}{\kappa ^2}, \end{equation} \tag{ 18 }$

$\begin{equation} {\tilde{\tau }}_{m} \equiv {\tau }_{m}\left(\frac{{\Sigma }_{m}}{\Sigma }\right)^{-1} = \frac{c_{ s}^2}{{\eta }_{m}} \int ^r_0 \frac{{\eta }_{m}}{c_{ s}^2}\frac{d\tau }{dr^{\prime }} dr^{\prime }. \end{equation} \tag{ 19 }$

$\mathscr{R}$ has the same units as the inverse of vortensity, but is a quantity that depends on radiation pressure. ${\tilde{\tau }}_{m}$ is the ratio between the perturbed optical depth τ_m and the relative surface density perturbation Σ_m/Σ. The local disk is unstable if a solution for k_r exists given ω² = 0, which denotes the line of neutral stability. Setting ω² = 0, the condition for $k_{ r}^2>0$ is

$\begin{equation} \beta e^{-\tau } \left(\frac{\kappa }{\Omega _{ k}}\right)^{-2} \left(\frac{d\tau }{d\ln {r}} + {\tilde{\tau }}_{m}\frac{d\ln {[r\mathscr{R}]}}{d\ln {r}} \right) > 1. \end{equation} \tag{ 20 }$

It is important to note that κ contains dependences on both radiation and gas pressure. In the interest of specifically studying IRI, we consider the case when the rotation curve is solely modified by radiation pressure. Then κ can be expressed as

$\begin{equation} \left(\frac{\kappa }{\Omega _{ k}}\right)^{2} = 1 - \beta e^{-\tau }\frac{d\ln {[r\beta ]}}{d\ln {r}} + \beta e^{-\tau } \frac{d\tau }{d\ln {r}}. \end{equation} \tag{ 21 }$

Plugging Equation (21) into Equation (17), the condition for instability becomes

$\begin{equation} q_{\beta } \equiv \beta e^{-\tau } \left(\frac{d\ln {[r\beta ]}}{d\ln {r}} + {\tilde{\tau }}_{m} \frac{d\ln {[r\mathscr{R}]}}{d\ln {r}} \right) > 1. \end{equation} \tag{ 22 }$

To complete our derivation, we need to evaluate ${\tilde{\tau }}_{m}$ . We begin by integrating Equation (19) by parts:

$\begin{equation} {\tilde{\tau }}_{m} = \tau - ik_{ r} \frac{c_{ s}^2}{{\eta }_{m}} \int ^r_0 \frac{{\eta }_{m}}{c_{ s}^2} \tau dr^{\prime }. \end{equation} \tag{ 23 }$

If in the disk there exists a "transition region" where the disk sharply transitions from being radially transparent to opaque, then one can show that inside this region, the second term on the right-hand side of Equation (23) has a magnitude of the order of k_rΔr, where Δr is the width of the transition region. This allows us to approximate ${\tilde{\tau }}_{m} \sim \tau$ in the limit 1/Δr ≫ k_r. Moreover, even when 1/Δr ∼ k_r, we expect ${\tilde{\tau }}_{m} \sim \tau$ to remain accurate to within the order of unity. In Section 5.3, we evaluate ${\tilde{\tau }}_{m}$ explicitly and find out to what extend this holds true.

While Equation (20) is the more general form, Equation (22) does reveal surprising behavior: it contains no explicit dependence on dτ/dr, as it is completely canceled by the stabilizing effect of κ². Replacing it is a term containing dβ/dr, whose effect is to lower κ² to the point of triggering a form of irradiation-induced Rayleigh instability. While it does contribute to the instability of the disk, we do not consider it the true trigger of IRI. Rather, we focus on the second term inside the bracket. First, it implies that a disk is unstable to IRI if it has a positive gradient in $\mathscr{R}$ , which can be created by a gradient in Σ and/or β. Second, this gradient must be located where ${\tilde{\tau }}_{m}e^{-\tau }\sim \tau e^{-\tau }$ is reasonably large, which is precisely the transition region. This is consistent with our picture that IRI is driven by shadowing. Because of the uncertainty in ${\tilde{\tau }}_{m}$ , as well as the other assumptions stated in the beginning of this section, Equations (20) and (22) should be taken as order-of-magnitude guidelines rather than rigid conditions.

2.2. Corotating Modes

If a linear mode exists, its corotation radius can be found by solving Equation (16) for Ω_m = 0. Similar to Section 2.1, we apply the WKB approximation, and then the real part of Equation (16) evaluated at the corotation radius can be rewritten as

$\begin{equation} {\Omega }_{m} = 0 = 2m\Omega \frac{\frac{h^2}{r^2} \frac{d\ln {\mathscr{F}}}{d\ln {r}} - \beta e^{-\tau }{\tilde{\tau }}_{m}}{ \frac{\kappa ^2}{\Omega _k^2} + |k|^2 h^2 - \beta e^{-\tau } \left(\frac{d\tau }{d\ln {r}} + {\tilde{\tau }}_{m} \frac{d\ln {[r\mathscr{R}]}}{d\ln {r}} \right) }, \end{equation} \tag{ 24 }$

where $|k|^2 = k_r^2 + m^2 / r^2$ , and $\mathscr{F} \equiv \Sigma \Omega / \kappa ^2$ is a quantity inversely proportional to the vortensity of the disk. The corotation radius is therefore located where the following condition is satisfied:

$\begin{equation} \frac{d\ln {\mathscr{F}}}{d\ln {r}} = \left(\frac{h}{r}\right)^{-2}\beta e^{-\tau }{\tilde{\tau }}_{m}. \end{equation} \tag{ 25 }$

For barotropic flow and β = 0, this condition becomes identical to that described in Section 2.2 of Lovelace et al. (1999) for RWI. The usefulness of Equation (25) is limited because without a full solution, the exact value of ${\tilde{\tau }}_{m}$ is unknown. However, allowing that ${\tilde{\tau }}_{m} \sim \tau$ , it does provide an insight: since the right-hand side of Equation (25) is always positive, if $\mathscr{F}$ contains a local maximum, the corotation radius will always be located at a lower orbit than where this maximum is. In our disk model described in the following section, $\mathscr{F}$ does contain a local maximum within the transition region, so we expect the corotation radius to be smaller for disks with a larger value of (h/r)⁻²β. This prediction is tested in Section 5.1.

3. DISK MODEL

For simplicity, we do not consider any spatial variation in the composition of the disk; therefore β and κ_opa are constants. With this simplification, Equation (22) says that the disk is most unstable if $\mathscr{R}$ has a large positive gradient near τ = 1. We create this condition with a disk that contains a sharp inner edge. At this edge, Σ increases by orders of magnitude across a small radial range, while τ rises from a small value to above unity. Our prescription for such a disk is

$\begin{equation} \Sigma (r) = \frac{1}{2}(\Sigma _{d}+\Sigma _{ c})+\frac{1}{2}{\rm erf}\left(\frac{r-r_0}{\sqrt{2{\Delta r}^2}}\right)(\Sigma _{d}-\Sigma _{ c}), \end{equation} \tag{ 26 }$

where Σ_d is the surface density of the disk, Σ_c is the surface density inside the cavity, r₀ is the radius at which the inner edge is located, and Δr is the width of this edge. We set Σ_d = 1 and Σ_c = 0.001 for a density contrast of 10³. We also set r₀ = 1 and GM = 1 so that the dynamical time t_dyn at the edge is $\Omega _{ k}^{-1}=1$ . For the sharpness of the edge, we set Δr = 0.05. The motivation for this choice is that Δr is unlikely to be shorter than h, which for protoplanetary disks has a typical value of 0.05r. κ_opa is chosen such that τ(r₀) = 1. If we move this τ = 1 point to a much smaller/larger radius, the disk edge will be become optically thick/thin, and thus one would expect the instability to weaken or even disappear. Figure 2 plots both the Σ and τ profile. To complete the equation set, we adopt an isothermal equation of state so that c_s is a constant.

**Figure 2.** Black solid line plotting the surface density profile described by Equation (26) and red dashed line plotting the optical depth profile.
Download figure:
Standard image High-resolution image

This leaves two free parameters in our model: β and c_s. We perform a parameter study over the range β = {0, 0.3} and c_s = {0.02, 0.06}. Note that for h(r₀) ≳ Δr, corresponding to c_s ≳ 0.05, the disk edge may become hydrodynamically unstable. We deliberately include this limit in our parameter space both as a sanity check and to investigate how IRI can be differentiated from other forms of instabilities.

4. TWO INDEPENDENT APPROACHES

For our given disk model, we aim to find out for the IRI (1) how the modal growth rate varies as a function of β and c_s, and (2) what are the properties of its nonlinear phase. Two independent approaches are used: a numerical method using hydrodynamical simulations and a semi-analytic method that solves the linearized problem (Equation (16)). These two methods not only serve as verifications for each other, but are also complementary since a full simulation gives us an insight into the nonlinear phase, while the semi-analytic method is not subjected to limitations such as resolution and numerical noise.

4.1. Hydrodynamical Simulation

We numerically simulate the 2D disk described in Section 3. The code we use is the Lagrangian, dimensionally split, shock-capturing hydrodynamics code PEnGUIn (Piecewise Parabolic Hydro-code Enhanced with Graphics Processing Unit Implementation), which has been previously used to simulate disk gaps opened by massive planets (Fung et al. 2014). It uses the piecewise parabolic method (PPM; Colella & Woodward 1984), and its main solver is modeled after VH-1 (Blondin & Lufkin 1993), with a few of the same modifications as described in Fung et al. (2014). It solves Equations (2) and (3), and contains an additional module to compute τ using piecewise parabolic interpolation to match the order of the PPM.

Our simulations have a domain spanning 0.5–2.0 in radial (in units where the disk edge is located at r₀ = 1) and the full 0–2π in azimuth. Moving the inner boundary to 0.7 or the outer boundary to 1.5 has a negligible effect on the growth of linear modes. We opt for a larger domain to accommodate the more violent nonlinear evolution.

The resolution is 1024 (r) by 3072 (ϕ). Azimuthal grid spacing is uniform everywhere, but radial grid spacing is uniform only between 0.5 and 1.3; from 1.3 to 2.0 it is logarithmic. This takes advantage of PEnGUIn's ability to utilize non-uniform grids to enhance the resolution around the disk edge. The resulting grid size at r₀ is about 0.001 (r) by 0.002 (ϕ). This gives at least 10 cells per h for even the smallest h we consider. Figure 3 shows how our simulations converge with resolution.

**Figure 3.** Growth rates of azimuthal modes with (β, *c_s*) = (0.2, 0.02). At 1024 (r) by 3072 (ϕ), the growth rates extracted from simulation match those found by the semi-analytic method to ∼1%. For this particular case, the fastest growing mode is m = 18, with a growth rate of ${\rm Im}(\omega)=4.0\times 10^{-1} t_{\rm dyn}^{-1}$ . See Section 5.1 for further discussions on how these results vary with β and c_s.
Download figure:
Standard image High-resolution image

**Figure 3.** Growth rates of azimuthal modes with (β, *c_s*) = (0.2, 0.02). At 1024 (r) by 3072 (ϕ), the growth rates extracted from simulation match those found by the semi-analytic method to ∼1%. For this particular case, the fastest growing mode is m = 18, with a growth rate of ${\rm Im}(\omega)=4.0\times 10^{-1} t_{\rm dyn}^{-1}$ . See Section 5.1 for further discussions on how these results vary with β and c_s.
Download figure:
Standard image High-resolution image

In each simulation, we extract the amplitudes of azimuthal modes as functions of time, resolving up to m = 50:

$\begin{equation} {A}_{m}(t) = \frac{1}{2\pi }\left|\int ^{2\pi }_{0}\int ^{1.1}_{1.0}\Sigma (t) e^{im\phi } dr d\phi \right|, \end{equation} \tag{ 27 }$

where we have chosen to integrate over the radial range r = {1.0, 1.1}. Instantaneous values of A_m are not the focus; rather, we seek a distinct period of exponential growth where we can measure its growth rate, i.e., the imaginary part of ω. Figure 4 shows one example of how modal growth behaves in these simulations. For a disk of a given set of parameters, the highest growth rate characterizes its timescale for instability.

**Figure 4.** Temporal evolution of *A_m* (see Equation (27)) with (β, *c_s*) = (0.05, 0.05). A well-defined exponentially growing phase can be seen around t = 200 ∼ 300. Beyond t = 300, the modes begin to exhibit higher-order coupling.
Download figure:
Standard image High-resolution image

We use a boundary condition fixed to the initial values described by Equations (5) and (26), with zero radial velocity. To reduce noise in A_m, we also include wave-killing zones in r = {0.5, 0.6} for the inner boundary and r = {1.6, 2.0} for the outer. Within these zones, we include an artificial damping term:

$\begin{equation} \frac{\partial {X}}{\partial t} = \left(X(t=0)-X\right)\frac{2c_{ s}|r-r_{\rm kill}|}{d_{\rm kill}^2}, \end{equation} \tag{ 28 }$

where X includes all disk variables Σ, P, and $\bf v$ ; r_kill is the starting radius of the wave-killing zone, which is 0.6 for the inner boundary and 1.6 for the outer; and d_kill is the width of these zones, which equals 0.1 for the inner boundary and 0.4 for the outer. In the end we are able to resolve A_m as small as 10⁻¹⁰, such as shown in Figure 4.

Simulations are terminated soon after the instability becomes fully nonlinear: up to 100 orbits, or 628 t_dyn. For very slowly growing modes, numerical noise severely hampers the precision of growth rate measurements. Consequently, this method is only capable of measuring growth rates ${\gtrsim } 0.01 t_{\rm dyn}^{-1}$ . The computational time for PEnGUIn is about 12 minutes per orbit on a single GTX-Titan graphics card.

4.2. Semi-analytic Method

Equation (16) constitutes an eigenvalue problem, where η_m is the eigenfunction and ω is the eigenvalue. To solve this problem, we develop a code that directly integrates the differential equations, iterates for the correct boundary conditions, and optimizes to find the eigenvalues. The complexity of this code is mainly to overcome the difficulty imposed by the integral in Equation (16), which effectively raises the order of the differential equation. The details are documented in the Appendix.

Despite the fact that it solves the linearized equations, our semi-analytic method in fact requires a much longer computational time than simulations using PEnGUIn. Due to limited resources, initially we only apply it to five sets of parameters. Table 1 contains a list of these sets. One major advantage of this method is that it does not have a limit to how slow of a growth rate can be detected, so we also apply it to all cases where simulations do not detect any modal growth. Among them, we find a positive growth rate for one case.³ It is also listed in Table 1, making a total of six sets of parameters.

Table 1. Semi-analytic Results^a

β	c_s	Im(ω)	m	r_cor^b
β	c_s	$(t_{\rm dyn}^{-1})$	m	(r₀)
0	0.06	7.7 × 10⁻²	4	1.046
0	0.05	2.5 × 10⁻³	1	1.052
0.1	0.05	7.2 × 10⁻²	6	0.979
0.15	0.04	1.6 × 10⁻²	7	0.956
0.15	0.03	1.4 × 10⁻¹	8	0.943
0.2	0.02	4.0 × 10⁻¹	18	0.938

Notes. ^aWe only report the properties of the fastest growing mode. ^br_cor denotes the corotation radius.

Download table as: ASCII Typeset image

5. RESULTS

The sets of parameters we consider are β = {0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3} and c_s = {0.02, 0.03, 0.04, 0.05, 0.06}. All 30 combinations of these values are simulated, but only a select few are solved with our semi-analytic method (see Table 1). The growth rates found by our two independent approaches agree to ∼1%. Figure 3 gives one example of this agreement. Also, the shapes of the modes extracted from simulations are nearly identical to the ones solved semi-analytically. Comparing Figures 5 and 6, results from the two methods are only distinguishable near the outer boundary, where the simulated ones show some artificial damping due to the wave-killing zones imposed. Because of the excellent agreement, we are able to combine the results of the two approaches to give a detailed picture for the IRI linear modes, complemented by the nonlinear evolution provided by simulations.

**Figure 5.** Fastest growing modes extracted from simulations through Fourier decomposition. Color shows the surface density normalized to the peak of each mode. On the left is an m = 18 mode from (β, *c_s*) = (0.2, 0.02); in the middle is m = 6 from (β, *c_s*) = (0.1, 0.05); and on the right is m = 4 from (β, *c_s*) = (0, 0.06).
Download figure:
Standard image High-resolution image

**Figure 6.** Fastest growing modes directly computed using our semi-analytic method for the same parameters listed in Figure 5.
Download figure:
Standard image High-resolution image

5.1. Linear Modes

We find clear growth of asymmetric modes for all cases with β larger than a certain threshold value that is weakly dependent on c_s across our parameter space. For most of our chosen c_s values, modal growth is only detected when β ⩾ 0.1, except for c_s ∼ 0.02, where this threshold rises to β ⩾ 0.15. From a simple perspective, we expect the disk to be more unstable for a larger β and smaller c_s, because β measures the strength of radiation pressure while c_s is a source of resistance to external forcing. In general, we do find the growth rate to increase with β and decrease with c_s, but with obvious exceptions.

In Figure 7, we divide our parameter space into three regions: regions I and II where modal growth is driven by radiation pressure, and region III where it is mainly driven by hydrodynamical effects. In regions I and II, growth rate scales roughly linearly with β for any given c_s, a trend that can be more easily seen in Figure 8. This is consistent with our expectations.

**Figure 8.** Growth rate of the fastest growing mode as a function of β. The (β, *c_s*) = (0, 0.05) point is disconnected because no modal growth is detected at (β, *c_s*) = (0.05, 0.05).
Download figure:
Standard image High-resolution image

Figure 9, however, reveals a more complicated aspect of IRI. Disregarding the β = 0 data points that belong to region III, the growth rate is very close to constant over the range 0.04 ⩽ c_s ⩽ 0.06 for any given β. Once c_s goes below 0.04 it shows different trends depending on the value of β: the growth rate increases as c_s decreases for β ⩾ 0.2, but for a smaller β the trend flattens or even begins to drop. This complex behavior may relate to how sound waves and IRI modes couple. While sound waves have a length scale h, IRI modes are mainly restricted by the sharpness of the transition region, which has a length scale Δr. Our results suggest that the coupling is weak when h ∼ Δr, and becomes much stronger as h decreases, allowing the transition region to accommodate the full wavelength of the longest wave.

Region III is where radiation pressure becomes a smaller effect than gas pressure. For c_s ⩾ 0.05, or equivalently, h(r₀) ⩾ Δr, we detect modal growth even in the absence of any radiation pressure. In fact, when β = 0, our disk model is similar to the "homentropic step jump" model used by Li et al. (2000) to study RWI (see their Figure 2). Their Figure 11 shows that for a pressure jump with a width Δr = 0.05, RWI modes will develop if c_s ≳ 0.06.⁴ Our results are consistent with their findings.

The division between IRI and RWI is clear to us because the two mechanisms appear to destructively interfere with each other. For c_s = 0.06, there is a clear drop in growth rate from β = 0 to β = 0.05 before it rises again (see Figure 8). Similarly, for c_s = 0.05, we do not detect any modal growth at β = 0.05 even though it is detected at both β = 0 and β = 0.1. One clue to this behavior is that we find c_s and β to have opposing effects on the epicyclic frequency κ. In Figure 10, we see that gas pressure lowers κ near r = r₀, while β raises it. Two effects roughly cancel when β = 0.05. It is unclear whether this is coincidental or not.

**Figure 10.** κ vs. r for three sets of parameters. The black dotted curve shows κ modified by gas pressure only. As β increases, the local minimum near r = 1 is flattened (red solid curve) and reversed (blue dashed curve).
Download figure:
Standard image High-resolution image

This dividing line may not remain at β = 0.05 for a different value of Δr or c_s. If we create a sharper edge by reducing Δr, both IRI and RWI are expected to be enhanced and it is unclear to us whether this dividing line will move to a higher or lower β. Additionally, a high c_s can push κ² below zero and trigger Rayleigh instability, which further complicates the matter. For our disk model, this limit is at c_s ∼ 0.07. Since the focus of this paper is to characterize IRI, we defer the thorough investigation on the interactions between IRI and other forms of instability to a future study.

Other than the growth rate, we also find other general trends about the linear modes. For an increasing β or decreasing c_s, the azimuthal mode number m of the fastest growing mode increases. The dependence on β is particularly pronounced. In the most extreme case, the fastest growing mode is m = 47 when (β, c_s) = (0.3, 0.02). In contrast, for all cases with β = 0.1, m = 5 ∼ 6 is the fastest growing mode. Figure 3 shows an intermediate case where the fastest growing mode is m = 18. See Table 1 for more examples. Similarly, we find that the radial extent of the fastest growing mode becomes more confined as β increases and c_s decreases, as can be seen in Figures 5 and 6. It is therefore empirically apparent that a higher β encourages the growth of a shorter wavelength mode. While our theory does not make any predictions about which mode grows the fastest, in hindsight this result is not surprising because the radial motion of an IRI mode must be driven by radiation pressure, so a stronger perturbing force should generate a faster radial motion, and therefore a higher frequency wave.

Another trend is that for an increasing β or decreasing c_s, the corotation radius decreases. This is in accordance with our prediction in Section 2.2. The dependence is weak but noticeable (see Table 1). Curiously, Figure 5 and 6 show that the peak location of each mode is relatively insensitive to variations in both β and c_s. Consequently, these peaks generally do not coincide with their corotation radii.

For the bulk of this work, we do not explicitly vary Δr as a free parameter. Since our choice of Δr = 0.05 is arbitrary, it is useful to find out to what extent our results would change for a different Δr. Setting c_s = 0.04 and Δr = 0.1, we find the threshold for modal growth becomes β ⩾ 0.25, roughly twice as large as when Δr = 0.05. This suggests that the threshold value scales roughly linearly with Δr for the parameter space we considered.

5.2. Nonlinear Evolution

Nonlinear evolution is what separates region I from region II⁵ in Figure 7. Figures 11 and 12 shows the simulation snapshots for the same two sets of parameters in the left and middle panels of Figures 5 and 6. The left panel, which belongs to region I with (β, c_s) = (0.2, 0.02), shows local regions of very high surface density, exceeding 10 times Σ_d, the initial surface density of the disk defined in Equation (26). This type of "clumping," which we define as a detection of Σ > 2Σ_d anywhere in the disk, is characteristic of region I. Note that clumping is not a necessary product of IRI, since region II is also driven by IRI. The transition from region II to I is rapid, in the sense that as we move toward the upper left corner of Figure 7, the highest local surface density quickly rises to a few tens of Σ_d.

**Figure 12.** Cartesian view of Figure 11.
Download figure:
Standard image High-resolution image

To compare the two regions in detail, we use the right panel of Figure 11 as a typical case for region II. It shows a significant widening of the edge by a factor of ∼ three. Vortices are formed along the edge and they create mild local enhancements in density. Their structure is complex as they typically launch two sets of spiral arms instead of one. The number of vortices is initially equal to the mode number of the fastest growing linear mode, but as they interact with each other, they can occasionally merge. It is unclear how many will remain in the long run since our simulations only last for 100 orbits at most.

In comparison to region II, the clumping in region I creates a much different, almost violent, nonlinear evolution. The clumps are very sharp features, with jumps over three orders of magnitude in density while their sizes are merely ∼0.1r₀. They are constantly formed and destroyed by disk shear over a dynamical timescale. The destroyed clumps form high density streams that are also visible in the left panel of Figure 11. Another consequence of this clumping is that by concentrating a large amount of matter in a small region, radiation is able to penetrate further into the disk and push the edge of the disk to a higher orbit. In the same figure, one can see that the edge of the disk is shifted to ∼1.2r₀. We speculate that the difference between regions I and II is due to the influence of gas pressure. Higher gas pressure results in vortex formation more similar to the purely hydrodynamical RWI, and as gas pressure becomes weaker compared to radiation, sharper features are created.

5.3. ${\tilde{\tau }}_{m}$ and the Instability Criterion Revisited

In our derivation for the instability criterion in Section 2.1, we propose the crude assumption ${\tilde{\tau }}_{m}\sim \tau$ . Using our semi-analytic method to obtain solutions for η_m, we are able to evaluate ${\tilde{\tau }}_{m}$ explicitly. For all IRI modes we have solved semi-analytically, we find ${\tilde{\tau }}_{m}/\tau >1$ within the region r = {r₀ − Δr, r₀}, but it is never an order of magnitude above unity. For example, Figure 13 plots ${\tilde{\tau }}_{m}/\tau$ for the m = 18 mode with (β, c_s) = (0.2, 0.02). It shows that within 0.93 ⩽ r ⩽ 1.02. The real part of ${\tilde{\tau }}_{m}/\tau$ is within 1 to 6, while the imaginary part is close to vanishing. Using this empirical result we can rewrite Equation (22) as

$\begin{equation} q_{\beta } \approx \beta e^{-\tau } \left(\frac{d\ln {[r\beta ]}}{d\ln {r}} + f\tau \frac{d\ln {[r\mathscr{R}]}}{d\ln {r}} \right), \end{equation} \tag{ 29 }$

where we substitute ${\tilde{\tau }}_{m}$ for fτ, and f > 1 is a number of the order of unity. Choosing f = 3, Figure 14 plots q_β for a few different β. This choice of f puts the threshold for instability (q_β > 1) at around β = 0.1, similar to our empirical results.

**Figure 13.** ${\tilde{\tau }}_{m}/\tau$ for the m = 18 mode with (β, *c_s*) = (0.2, 0.02). Inside the transition region, r ≈ {0.95, 1.0}, the approximation ${\tilde{\tau }}_{m}\sim \tau$ is accurate to within the order of unity.
Download figure:
Standard image High-resolution image

**Figure 13.** ${\tilde{\tau }}_{m}/\tau$ for the m = 18 mode with (β, *c_s*) = (0.2, 0.02). Inside the transition region, r ≈ {0.95, 1.0}, the approximation ${\tilde{\tau }}_{m}\sim \tau$ is accurate to within the order of unity.
Download figure:
Standard image High-resolution image

**Figure 14.** q_β from Equation (29) for different values of β. We choose f = 3 to best match our empirical results.
Download figure:
Standard image High-resolution image

q_β becomes negative as soon as τ > 1 (r > r₀ = 1 for our disk model), because the exponential factor in $\mathscr{R}$ quickly drives its gradient negative. This implies that the instability must originate from the τ < 1 region. Along the same line, we find the corotation radii of IRI modes to be within r₀, as shown in Table 1.

6. CONCLUSIONS AND DISCUSSIONS

We demonstrated that IRI can operate at an inner disk edge where there is a transition from being radially transparent to opaque. A local criterion for axisymmetric instability was derived (Equation (22)). For our given disk model, we computed the linear modal growth rates for β varying from 0 to 0.3, and c_s from 0.02 to 0.06. We found growth rates ranging from 10⁻² to $10^0 t_{\rm dyn}^{-1}$ (Figure 7). The fastest rates were found for the largest β and smallest c_s. We empirically determined that the threshold for IRI is β ∼ 0.1 when Δr = 0.05, with a weak dependence on c_s. For a wider edge, Δr = 0.1, this threshold rises to β ∼ 0.25. We note that this implies the threshold can be lowered by reducing Δr; however, at the same time c_s must also be lowered for IRI to dominate over other forms of instability that may be triggered by the sharpness of the edge, such as RWI and Rayleigh instability. We employed two independent approaches to obtain the growth rates of the linear modes: simulating the disks numerically using PEnGUIn, and solving the linearized equations semi-analytically. Their excellent agreement lends confidence in our results. Moreover, we discovered a parameter space, labeled region I in Figure 7, where "clumping" occurs. There one can find over 10 times the local surface density enhancements in the nonlinear evolution of IRI.

6.1. Connection to Physical Disks

Our disk model is inspired by transitional disks (e.g., Calvet et al. 2005; Espaillat et al. 2007, 2008; Andrews et al. 2011). The inner edges of these disks are currently unresolved by observation, but theoretical work has shown that the sharpness of disk edges created by X-ray photoevaporation (e.g., Owen et al. 2010) is similar to that described by our Equation (26) with Δr = 0.05 (compare our Figure 2 to Figure 2 of Owen et al. 2013). If a transitional disk undergoes IRI, the asymmetric structure at the inner edge will create an azimuthal variation in shadowing. Flaherty & Muzerolle (2010) showed that this can lead to a significant variation in disk emission. Indeed, some variability in the infrared emission of transitional disks has been reported by Muzerolle et al. (2009), Flaherty et al. (2011), and Espaillat et al. (2011).

On the other hand, IRI is by no means limited to circumstellar disks. AGN accretion disks, for example, can be subjected to IRI if there are any sharp jumps in density and/or opacity, such as the inner edges of the broad-line regions. IRI can potentially generate the stochastic asymmetry, which is used to explain the variability in the double-peaked Balmer emission lines in radio-loud AGNs (Flohic & Eracleous 2008). We note that the dynamics in AGN accretion disks are considerably more complicated since they do not have a point-like light source.

6.2. Implications of "Clumping"

The "clumping" found in a part of our parameter space (Figure 7) opens new possibilities for IRI. For instance, very high density regions in protoplanetary disks may be favorable environments for the formation of planetary cores. The density of individual clumps may even become high enough to trigger gravitational instability at the inner edges of massive disks. One should be cautious to interpret the enhancement factors reported as realistic, however, since it is only one disk model that we have studied.

The clumping also leads to a possibility of preventing inward dust migration. Dominik & Dullemond (2011) demonstrated that while radiation pressure can initially push dust outward and form a dust wall, the wall eventually succumbs to the global accretion flow and migrates inward. If this wall becomes unstable due to IRI, clumping can occur, effectively creating "leakage" within the wall, allowing radiation to push dust further back. The true behavior of these dust walls is important to understand disks where inner clearings have been observed, such as transitional disks. Dynamical interactions between radiation, dust, and gas must be considered for this kind of study.

6.3. Outlook

There are three main aspects of our model that we feel would benefit greatly from a more realistic treatment. First, our model ignores the vertical dimension. A notable difference from two dimensions to three dimensions is that the location of the inner edge of a disk, defined as the τ = 1 point, would become a function of height, spreading over a distance of ∼h. One possible consequence is that IRI would generate a vertical circulation at the inner edge, which would dilute the opacity in the midplane and allow radiation pressure to penetrate further into the disk. Additionally, in a flared disk, radiation pressure is exerted on the photosphere of the entire disk rather than just the inner edge. On the other hand, because of dust settling, we expect the value of β in the photosphere to be smaller than the midplane, making it even more difficult to reach the β ≈ 0.1 threshold. Nonetheless, for disks around exceptionally luminous stars, IRI can potentially operate at all radii.

Second, we assume a perfect coupling between gas and dust. In a more realistic approach, dust should be allowed to migrate with respect to gas. One expects dust to gather near the initial τ = 1 point, because where it is optically thin, dust migrates outward due to the effect of radiation pressure, and in the optically thick disk, dust migrates inward due to gas drag. This behavior of dust is described in Section 3 of Takeuchi & Artymowicz (2001). The buildup of a dust wall is almost certain to trigger IRI due to its large β gradient.

Lastly, we lack a realistic treatment for radiative transfer. As the disk crosses from being radially transparent to opaque, the midplane of the disk also transitions from being heated directly by irradiation, to passively by the irradiated atmosphere. Consequently, the midplane temperature should be decreasing across the disk edge. This is not captured by our globally isothermal assumption. Additionally, the clumps we find in some of our nonlinear results are sufficiently dense that they are optically thick. With our isothermal treatment, they remain the same temperature as their surroundings, while in truth these clumps should be capable of shielding themselves from irradiation and creating a non-trivial internal temperature structure. Whether this is an effect that aids or inhibits their formation and survival requires future investigation.

We thank Yanqin Wu, Chris Matzner, and Eugene Chiang for helpful feedback that substantially improved this manuscript. We also thank James Owen for insightful discussions. J.F. owes his gratitude to the Queen Elizabeth II Graduate Scholarship in Science and Technology. We gratefully acknowledge support from the Discovery Grant by the Natural Sciences and Engineering Research Council of Canada.

APPENDIX: NUMERICAL METHOD FOR SOLVING THE LINEARIZED EQUATIONS

In this Appendix, we document our method for solving Equation (16) numerically. To begin, note that Equation (16) can only be numerically integrated in the direction of increasing r because of the integral in the fourth term. In principle, it is possible to simply do this integration and find the value of ω that best matches the desired outer boundary condition. This is impractical, however, because any slight error in ω leads to a diverging behavior of η_m at the outer boundary. A better method is to integrate Equation (16) simultaneously from the inner boundary outward, and the outer boundary inward, and find the ω that results in a match of the two functions at some intermediate radius r_mid. To accomplish this, we first define

$\begin{equation} {y}_{m} \equiv \int ^r_0 \frac{{\eta }_{m}}{c_{ s}^2}\frac{d\tau }{dr^{\prime }}dr^{\prime }, \end{equation} \tag{ A1 }$

and then differentiate Equation (16) with respect to r:

$\begin{equation} \frac{\partial ^3{y}_{m}}{\partial r^3}+a^{\prime }(r)\frac{\partial ^2{y}_{m}}{\partial r^2}+b^{\prime }(r)\frac{\partial {y}_{m}}{\partial r}+c^{\prime }(r){y}_{m} = 0, \end{equation} \tag{ A2 }$

where

$\begin{eqnarray*} \nonumber a^{\prime } &\equiv a - 2\frac{d\ln {g}}{dr},\\ \nonumber b^{\prime } &\equiv b - a\frac{d\ln {g}}{dr} + 2\left(\frac{d\ln {g}}{dr}\right)^2 - \frac{1}{g}\frac{{\rm d}^2 g}{dr^2},\\ \nonumber c^{\prime } &\equiv c g,\\ \nonumber g &\equiv \frac{1}{c_{ s}^2}\frac{d\tau }{dr}. \end{eqnarray*}$

Thus, we can now numerically integrate Equation (A2) in both directions, and recover η_m from y_m. The boundary conditions can be approximated using the WKB method. The WKB form for y_m is

$\begin{equation} {y}_{m} = R(r)e^{i\int ^r_0 kdr^{\prime }}, \end{equation} \tag{ A3 }$

$\begin{equation} \frac{\partial {y}_{m}}{\partial r} \simeq ik{y}_{m}, \end{equation} \tag{ A4 }$

where R(r) is a slowly varying function and k is the complex wave number that satisfies |kr| ≫ 1. Substituting Equations (A3) and (A4) into Equation (A2), we get the following algebraic equation for k:

$\begin{equation} k^3 - i a^{\prime } k^2 - b^{\prime } k + ic^{\prime } = 0. \end{equation} \tag{ A5 }$

The three solutions of Equation (A5) correspond to the inward traveling (Re(k) < 0), outward traveling (Re(k) > 0), and a third solution that does not exist in the conventional WKB approximation. In fact, it has |kr| ≪ 1, effectively rendering y_m a constant, which violates the approximation of a tightly winding wave. To accommodate for this solution, we generalize Equation (A3) to allow for a constant offset:

$\begin{equation} {y}_{m} = R(r)e^{i\int ^r_0 kdr^{\prime }} + C. \end{equation} \tag{ A6 }$

Substituting this into Equation (A2), we obtain:

$\begin{equation} k^3 - i a^{\prime } k^2 - b^{\prime } k + ic^{\prime }\left(\frac{{y}_{m}}{{y}_{m}-C}\right) = 0. \end{equation} \tag{ A7 }$

In the optically thin and thick limits, c' becomes arbitrarily small, and since the |kr| ≪ 1 solution is already incorporated into the constant offset C, the last term can be dropped, giving back the usual quadratic form:

$\begin{equation} k^2 - i a^{\prime } k - b^{\prime } = 0, \end{equation} \tag{ A8 }$

which gives the expected incoming and outgoing solutions for tightly winding waves. We apply the radiative boundary condition, assuming no wave is entering the domain from the boundaries. The other unknowns remaining in Equation (A6) are R and C. For clarity, we will denote variables associated with the solution integrated from the inner boundary with the subscript "in," and those from the outer boundary with "out."

Recall that y_m is in fact the integral of the perturbation (Equation (A1)). At the inner boundary, this quantity is small since inward of the boundary there is only a traveling wave, so we set C_in = 0. We choose R_in = 1, while R_out and C_out are determined by the following iterative formulas:

$\begin{equation} R_{\rm out}^{i+1} = R_{\rm out}^{i} \frac{d^2 y_{m,\rm in}^{i}}{dr^2} \left(\frac{d^2 y_{m,\rm out}^{i}}{dr^2}\right)^{-1}, \end{equation} \tag{ A9 }$

$\begin{equation} C_{\rm out}^{i+1} = C_{\rm out}^{i} + y_{m,\rm in}^{i} - y_{m,\rm out}^{i}, \end{equation} \tag{ A10 }$

where i is the current iterative step, and the y_m and its derivatives are evaluated at r_mid. Convergence typically requires tens or even hundreds of iterations, which is the primary reason for the large amount of computational time required for this method. Last, we find the eigenvalue ω by minimizing the following function, evaluated at r_mid:

$\begin{eqnarray} &&f = \left(\frac{{\rm Re}\left(\frac{dy_{m,\rm out}}{dr}\right)-{\rm Re}\left(\frac{dy_{m,\rm in}}{dr}\right)}{{\rm max}\left[\left|{\rm Re}\left(\frac{dy_{m,\rm out}}{dr}\right)\right|, \left|{\rm Re}\left(\frac{dy_{m,\rm in}}{dr}\right)\right|\right]}\right)^2 \nonumber\\ &&+ \left(\frac{{\rm Im}\left(\frac{dy_{m,\rm out}}{dr}\right)-{\rm Im}\left(\frac{dy_{m,\rm in}}{dr}\right)}{{\rm max}\left[\left|{\rm Im}\left(\frac{dy_{m,\rm out}}{dr}\right)\right|, \left|{\rm Im}\left(\frac{dy_{m,\rm in}}{dr}\right)\right|\right]}\right)^2. \end{eqnarray} \tag{ A11 }$

We use an eighth-order Runge–Kutta method with adaptive step-size control for the numerical integration. We set r_mid = 1, the inner boundary at r_in = 0.3, and the outer at r_out = 4. Minimizing f is also very time consuming because we employ a random sampling method: first, we bracket the minimum within a range of likely values for the real and imaginary part of ω, then we randomly select ω within the chosen range, and narrow down the field by preferentially choosing values closer to where f is below a certain threshold. This time consuming method is ultimately superior to methods that involve descending along the gradient of f, because of the numerous local minima that exist.

IRRADIATION INSTABILITY AT THE INNER EDGES OF ACCRETION DISKS

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION