Irradiation Instability at the Inner Edges of Accretion Disks
DD raft version N ovember
16, 2018
Preprint typeset using L A TEX style emulateapj v. 5 / / IRRADIATION INSTABILITY AT THE INNER EDGES OF ACCRETION DISKS J effrey F ung ( 馮 澤 之 ) and P awel A rtymowicz Draft version November 16, 2018
ABSTRACTAn instability can potentially operate in highly irradiated disks where the disk sharply transitions from beingradially transparent to opaque (the ”transition region”). Such conditions may exist at the inner edges of transi-tional disks around T Tauri stars and accretion disks around AGNs. We derive the criterion for this instability,which we term the ”irradiation instability”, or IRI. We also present the linear growth rate as a function of β ,the ratio between radiation force and gravity, and c s , the sound speed of the disk, obtained using two methods:a semi-analytic analysis of the linearized equations and a numerical simulation using the GPU-accelerated hy-drodynamical code PEnGUIn . In particular, we find that IRI occurs at β ∼ . ∼ . r , and at higher β values if it is wider. This threshold value applies to c s ranging from 3% of theKeplerian orbital speed to 5%, and becomes higher if c s is lower. Furthermore, in the nonlinear evolution of theinstability, disks with a large β and small c s exhibit ”clumping”, extreme local surface density enhancementsthat can reach over ten times the initial disk surface density. Subject headings: accretion, accretion disks, hydrodynamics, instabilities, protoplanetary disks, radiation: dy-namics INTRODUCTION
Accretion disks are susceptible to a wide range of in-stabilities, including the magnetorotational instability (MRI)(Balbus and Hawley 1998), gravitational instability (Lin andPringle 1987; Gammie 2001), Papaloizou-Pringle instability(Papaloizou and 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 shear-ing (Urpin and Brandenburg 1998) are considered. These in-stabilities drive the evolution of disks by generating turbu-lence and creating complex, sometimes extreme, structures,such as the formation of planets in protoplanetary disks.Radiation pressure is a force generally present in all typesof accretion disks. Its e ff ect on accretion disks has been stud-ied in many di ff erent aspects, including driving disk winds inactive galactic nuclei (AGN) (e.g. Higginbottom et al. 2014),shaping particle size distributions in debris disks (Thebaultet al. 2014), and influencing the motions of the inner rims oftransitional disks (Chiang and Murray-Clay 2007; Dominikand Dullemond 2011). We demonstrate in this paper that ra-diation pressure can also cause a disk instability of its ownkind. In the following, we give a brief introduction to thisinstability before launching into the formal theoretical work.The strength of radiation pressure compared to gravity ismeasured by the number β : β = κ opa L π cGM , (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 speedof light and gravitational constant respectively. The key tothis instability is shadowing. As the front part of the diskgets pushed by radiation pressure, it also casts a shadow that [email protected] Department of Astronomy and Astrophysics, University of Toronto,50 St. George Street, Toronto, ON, Canada M5S 3H4 Department of Physical and Environmental Sciences, University ofToronto at Scarborough, 1265 Military Trail, Scarborough, ON, CanadaM1C 1A4 reduces the amount of radiation pressure on the material fur-ther out in the disk. In a 1D, radial picture, since radiationpressure always diminishes outward, the inner part of a diskalways feels a stronger push than the outer part, and the nete ff ect is therefore radial compression. In other words, any twoconcentric disk annuli would feel an attraction between themdue to the combined e ff ects of radiation pressure and shadow-ing.This 1D scenario does not easily extend to a 2D disk how-ever, because radiation pressure from a central source doesnot exert any azimuthal force. By the conservation of angu-lar momentum, when a disk element is perturbed radially, itwill oscillate at some epicyclic frequency. Figure 1 illustrateswhat e ff ect this oscillating element has on the disk. One cansee that disk material near the orbit of the perturbed elementwill experience a variation in shadowing along the azimuth.This variation creates a forcing that induces the unperturbedmaterial to follow the motion of the perturbed element. Theresult is a global collective motion that is capable of growingon its own. We term this phenomenon the ”irradiation insta-bility” (IRI), since it relies on irradiation by the central object.Because a larger β allows for a more rapid radial motion,its value is crucial for the survival of this collective motionagainst disk shear. In most systems, dust grains provide thelargest contribution to β . In circumstellar disks, micron-sizegrains can have β > β ∼ for A-type stars (e.g., Equation 10 of Kirchschlager and Wolf(2013)). Given that the gas-to-dust ratio is typically ∼ , β of a perfectly coupled gas + dust mixture may be of the orderof a few percent. Additionally, dust settling can enhance β inthe midplane by reducing the local gas-to-dust ratio, while theradial migration of dust results in size segregation (Thebaultet al. 2014), which can also enhance β at local radii. In othersystems where radiation pressure can drive significant massloss, such as AGN accretion disks, one would even expect β to exceed unity.This paper aims to provide a basic understanding of IRI, ofboth the conditions that trigger it, and its consequences. InSection 2, we present a theoretical foundation for IRI and de-rive its instability criterion. Section 3 contains our disk model. a r X i v : . [ a s t r o - ph . E P ] J u l Fung & Artymowicz F ig . 1.— Simple illustration describing IRI. The blue curve denotes the orbitof a perturbed disk element oscillating around its guiding center, denoted bythe dashed black line at r . The shaded area is where the disk sees the shadowcast by the perturbed element. The red arrows show the directions of radialforcing on the background disk relative to the average amount of radiationpressure received along r . These arrows are inward when they are in theshadow of the element, and outward when they are not. One can see that thebackground disk near r is forced in the direction of amplifying the initialperturbation. Section 4 describes our semi-analytic and numerical methods.Section 5 reports the modal growth rate as a function of β andthe sound speed c s of the disk, and gives a discussion on thenonlinear evolution of IRI. Section 6 concludes with an out-look for future work. THE LINEAR THEORY
We follow the method of Goldreich and Tremaine (1979),using similar notation, to derive the linear response of a 2Ddisk stirred by radiation pressure. We start with the continuityequation and the conservation of momentum: ∂ Σ ∂ t + ∇ · ( Σ v ) = , (2) ∂ v ∂ t + ( v · ∇ ) v = −∇ η − GM (1 − β e − τ ) r ˆ r , (3)where Σ is the surface density of the disk; v is the 2D velocityfield; η is the specific enthalpy such that ∇ η = ∇ P / Σ , where P is the vertically averaged gas pressure; and τ is the opticaldepth of the disk. We denote the Keplerian orbital frequencyas Ω k , and the sound speed c s is defined by the ideal gas law P = c Σ . τ depends on the density distribution by the follow-ing equation: τ = (cid:90) r κ opa ρ dr (cid:48) , (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 theazimuthal direction. Σ , η , and v can be separated into a background quantity(without any subscript) and a perturbed quantity (denoted bythe subscript ” m ”). We assume the background disk to be ax-isymmetric and in hydrostatic equilibrium so that v = (0 , r Ω ),where Ω = (cid:114) Ω k (1 − β e − τ ) + r d η dr , (5)and the components of the perturbed velocity are denoted as v m ≡ ( u , υ ). To simplify the notation, we also define the back-ground and perturbed radiation force as F and F m : F = r Ω k β e − τ , (6) F m = − F (cid:90) r η m c d τ dr (cid:48) dr (cid:48) . (7)For a small perturbation, it follows from Equations 2 and3 that the perturbed quantities are governed by the followinglinearized equations: ∂ Σ m ∂ t + ∇ ( Σ v m ) + ∇ ( Σ m v ) = , (8) ∂ v m ∂ t + ( v · ∇ ) v m + ( v m · ∇ ) v = − ∇ η m + F m ˆ r . (9)Without a loss of generality, we can assume a form of thesolution for the perturbed quantities Σ m , η m , u and υ : X m ( r , θ, t ) = X ( r ) e i ( m θ − ω t ) , (10)for some complex function X ( r ) and complex number ω ,while m is the azimuthal mode number. Substituting this forminto Equation 9, we find u = − iD (cid:34) m Ω r η m + Ω m (cid:32) ∂η m ∂ r − F m (cid:33)(cid:35) , (11) υ = D (cid:34) m Ω m r η m + (cid:32) Ω + r d Ω dr (cid:33) (cid:32) ∂η m ∂ r − F m (cid:33)(cid:35) . (12)The pattern rotation frequency Ω m and the coe ffi cient D aredefined as Ω m ≡ m Ω − ω , (13) D ≡ κ − Ω m , (14) κ = r d (cid:104) r Ω (cid:105) dr , (15)where κ is the epicyclic frequency of the unperturbed orbit.To solve for Σ m , or equivalently η m , we substitute Equation11 and 12 into Equation 8, giving ∂ η m ∂ r + a ( r ) ∂η m ∂ r + b ( r ) η m + c ( r ) (cid:90) r η m c d τ dr (cid:48) dr (cid:48) = , (16)where a ≡ ∂∂ r ln (cid:32) r Σ D (cid:33) , b ≡ m Ω r Ω m ∂∂ r ln (cid:32) ΣΩ D (cid:33) − m r + c (cid:32) F d τ dr − D (cid:33) , c ≡ F (cid:32) ∂∂ r ln (cid:32) r Σ FD (cid:33) − m Ω r Ω m (cid:33) . rradiation Instability at the Inner Edges of Accretion Disks 3We arrive at a second-order integro-di ff erential equation for η m . Instability Criterion
A local criterion for axisymmetric instability can be derivedfrom Equation 16. We apply the WKB approximation andwrite η m ∼ e i (cid:82) r k r dr (cid:48) , where k r (cid:29) 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 writtenas: ω = κ + k r c − Ω k β e − τ (cid:32) d τ d ln r + ˜ τ m d ln [ r R ]d ln r (cid:33) , (17)where R ≡ ΣΩ k β e − τ κ , (18)˜ τ m ≡ τ m (cid:32) Σ m Σ (cid:33) − = c η m (cid:90) r η m c d τ dr (cid:48) dr (cid:48) . (19) R has the same units as the inverse of vortensity, but is aquantity that depends on radiation pressure. ˜ τ m is the ratio be-tween the perturbed optical depth τ m and the relative surfacedensity perturbation Σ m / Σ . The local disk is unstable if a so-lution for k r exists given ω =
0, which denotes the line ofneutral stability. Setting ω =
0, the condition for k r > β e − τ (cid:32) κ Ω k (cid:33) − (cid:32) d τ d ln r + ˜ τ m d ln [ r R ]d ln r (cid:33) > . (20)It is important to note that κ contains dependences on both ra-diation and gas pressure. In the interest of specifically study-ing IRI, we consider the case when the rotation curve is solelymodified by radiation pressure. Then κ can be expressed as (cid:32) κ Ω k (cid:33) = − β e − τ d ln (cid:2) r β (cid:3) d ln r + β e − τ d τ d ln r . (21)Plugging Equation 21 into Equation 17, the condition for in-stability becomes q β ≡ β e − τ (cid:32) d ln (cid:2) r β (cid:3) d ln r + ˜ τ m d ln [ r R ]d ln r (cid:33) > . (22)To complete our derivation, we need to evaluate ˜ τ m . Webegin by integrating Equation 19 by parts˜ τ m = τ − ik r c η m (cid:90) r η m c τ dr (cid:48) . (23)If in the disk there exists a ”transition region” where thedisk sharply transitions from being radially transparent toopaque, then one can show that inside this region, the secondterm on the right-hand side of Equation 23 has a magnitude ofthe order k r ∆ r , where ∆ r is the width of the transition region.This allows us to approximate ˜ τ m ∼ τ in the limit ∆ r (cid:29) k r .Moreover, even when ∆ r ∼ k r , we expect ˜ τ m ∼ τ to remain ac-curate to within the order of unity. In Section 5.3 we evaluate˜ τ m explicitly and find out to what extend this holds true.While Equation 20 is the more general form, Equation 22does reveal surprising behavior: it contains no explicit depen-dence on d τ dr , as it is completely canceled by the stabilizinge ff ect of κ . Replacing it is a term containing d β dr , whose e ff ect is to lower κ to the point of triggering a form of irradiation-induced Rayleigh instability. While it does contribute to theinstability of the disk, we do not consider it the true trigger ofIRI. Rather, we focus on the second term inside the bracket.First, it implies that a disk is unstable to IRI if it has a positivegradient in R , which can be created by a gradient in Σ and / or β . Second, this gradient must be located where ˜ τ m e − τ ∼ τ e − τ is reasonably large, which is precisely the transition region.This is consistent with our picture that IRI is driven by shad-owing. Because of the uncertainty in ˜ τ m , as well as the otherassumptions stated in the beginning of this section, Equations20 and 22 should be taken as order-of-magnitude guidelinesrather than rigid conditions. Corotating Modes
If a linear mode exists, its corotation radius can be foundby solving Equation 16 for Ω m =
0. Similar to Section 2.1,we apply the WKB approximation, and then the real part ofEquation 16 evaluated at the corotation radius can be rewrittenas: Ω m = = m Ω h r d ln F d ln r − β e − τ ˜ τ m κ Ω k + | k | h − β e − τ (cid:16) d τ d ln r + ˜ τ m d ln [ r R ]d ln r (cid:17) , (24)where | k | = k r + m / r , and F ≡ ΣΩ /κ is a quantity in-versely proportional to the vortensity of the disk. The corota-tion radius is therefore located at where the following condi-tion is satisfied: d ln F d ln r = (cid:32) hr (cid:33) − β e − τ ˜ τ m . (25)For barotropic flow and β =
0, this condition becomes iden-tical to that described in Section 2.2 of Lovelace et al. (1999)for RWI. The usefulness of Equation 25 is limited becausewithout a full solution, the exact value of ˜ τ m is unknown.However, allowing that ˜ τ m ∼ τ , it does provide an insight:since the the right-hand side of Equation 25 is always pos-itive, if F contains a local maximum, the corotation radiuswill always be located at a lower orbit than where this maxi-mum is. In our disk model described in the following section, F does contain a local maximum within the transition region,so we expect the corotation radius to be smaller for disks witha larger value of (cid:16) hr (cid:17) − β . This prediction is tested in Section5.1. DISK MODEL
For simplicity we do not consider any spatial variation inthe composition of the disk, therefore β and κ opa are constants.With this simplification, Equation 22 says that the disk is mostunstable if R has a large positive gradient near τ =
1. Wecreate this condition with a disk that contains a sharp inneredge. At this edge, Σ increases by orders of magnitude acrossa small radial range, while τ rises from a small value to aboveunity. Our prescription for such a disk is: Σ ( r ) =
12 ( Σ d + Σ c ) +
12 erf r − r (cid:112) ∆ r ( Σ d − Σ c ) , (26)where Σ d is the surface density of the disk, Σ c is the surfacedensity inside the cavity , r is the radius at which the inneredge is located, and ∆ r is the width of this edge. We set Σ d = Σ c = .
001 for a density contrast of 10 . We also set r = Fung & Artymowicz r -3 -2 -1 Σ || τ Στ F ig . 2.— Black solid line plotting the surface density profile described byEquation 26 and red dashed line plotting the optical depth profile. GM = t dyn at the edge is Ω − k =
1. For the sharpness of the edge, we set ∆ r = .
05. Themotivation for this choice is that ∆ r is unlikely to be shorterthan h , which for protoplanetary disks has a typical value of0 . r . κ opa is chosen such that τ ( r ) =
1. If we move this τ = / larger radius, the disk edge willbe become optically thick / thin, and thus one would expect theinstability to weaken or even disappear. Figure 2 plots boththe Σ and τ profile. To complete the equation set, we adopt anisothermal equation of state so that c s is a constant.This leaves two free parameters in our model: β and c s . Weperform a parameter study over the range β = { , . } and c s = { . , . } . Note that for h ( r ) (cid:38) ∆ r , correspondingto c s (cid:38) .
05, the disk edge may become hydrodynamicallyunstable. We deliberately include this limit in our parameterspace both as a sanity check and to investigate how IRI can bedi ff erentiated from other forms of instabilities. 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. Twoindependent approaches are used: a numerical method us-ing hydrodynamical simulations and a semi-analytic methodthat solves the linearized problem (Equation 16). These twomethods not only serve as verifications for each other, but arealso complementary since a full simulation gives us an insightinto the nonlinear phase, while the semi-analytic method isnot subjected to limitations such as resolution and numericalnoise. Hydrodynamical Simulation
We numerically simulate the 2D disk described in Section3. The code we use is the Lagrangian, dimensionally-split,shock-capturing hydrodynamics code
PEnGUIn ( P iecewiseParabolic Hydro-code En hanced with G raphics Processing U nit I mplementatio n ), which has been previously used tosimulate disk gaps opened by massive planets (Fung et al.2014). It uses the piecewise parabolic method (PPM; Colellaand Woodward 1984), and its main solver is modeled afterVH-1 (Blondin and Lufkin 1993), with a few of the samemodifications as described in Fung et al. (2014). It solvesEquation 2 and 3, and contains an additional module to com-pute τ using piecewise parabolic interpolation to match theorder of PPM.
10 15 20 25 30 35m0.20.250.30.350.4 I m ( ω ) [ t dyn - ] semi-analytic solutionsimulation (1024x3072)simulation (512x1536)simulation (256x768)simulation (128x384) F ig . 3.— Growth rates of azimuthal modes with ( β, c s ) = (0 . , . φ ), the growth rates extracted from simulation match those foundby the semi-analytic method to ∼ m =
18, with a growth rate of Im( ω ) = . × − t − . SeeSection 5.1 for further discussions on how these results vary with β and c s . Our simulations have a domain spanning 0.5 to 2.0 in radial(in units where the disk edge is located at r =
1) and the full0 to 2 π in azimuth. Moving the inner boundary to 0.7 or theouter boundary to 1.5 has a negligible e ff ect on the growth oflinear modes. We opt for a larger domain to accommodate themore violent nonlinear evolution.The resolution is 1024 (r) by 3072 ( φ ). Azimuthal gridspacing is uniform everywhere, but radial grid spacing is uni-form only between 0.5 and 1.3; from 1.3 to 2.0 it is loga-rithmic. This takes advantage of PEnGUIn ’s ability to utilizenon-uniform grids to enhance the resolution around the diskedge. 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 convergewith resolution.In each simulation, we extract the amplitudes of azimuthalmodes as functions of time, resolving up to m = A m ( t ) = π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) π (cid:90) . . Σ ( t ) e im φ d r d φ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (27)where we have chosen to integrate over the radial range r = { . , . } . Instantaneous values of A m are not the focus; rather,we seek a distinct period of exponential growth where we canmeasure its growth rate, i.e., the imaginary part of ω . Figure4 shows one example of how modal growth behaves in thesesimulations. For a disk of a given set of parameters, the high-est growth rate characterizes its timescale for instability.We use a boundary condition fixed to the initial values de-scribed by Equations 5 and 26, with zero radial velocity. Toreduce noise in A m , we also include wave-killing zones in r = { . , . } for the inner boundary and r = { . , . } forthe outer. Within these zones, we include an artificial damp-ing term: ∂ X ∂ t = ( X ( t = − X ) 2 c s | r − r kill | d , (28)where X includes all disk variables Σ , P , and v ; r kill is thestarting radius of the wave-killing zone, which is 0.6 for theinner boundary and 1.6 for the outer; and d kill is the width ofthese zones, which equals 0.1 for the inner boundary and 0.4for the outer. In the end we are able to resolve A m as small as10 − , such as shown in Figure 4.rradiation Instability at the Inner Edges of Accretion Disks 5
100 150 200 250 300 350 400 450time [ t dyn ]10 -12 -10 -8 -6 -4 -2 A m m=2m=3m=4m=5m=6m=7 F ig . 4.— Temporal evolution of A m (see Equation 27) with ( β, c s ) = (0 . , . t = ∼ t =
300 the modes begin to exhibit higher-ordercoupling. TABLE 1S emi - analytic R esults a β c s Im( ω ) ( t − ) m r corb ( r )0 0.06 7 . × − . × − . × − . × − . × − . × −
18 0.938 a We only report the properties of the fastestgrowing mode. b r cor denotes the corotation radius. Simulations are terminated soon after the instability be-comes fully nonlinear: up to 100 orbits, or 628 t dyn . For veryslowly growing modes, numerical noise severely hampers theprecision of growth rate measurements. Consequently, thismethod is only capable of measuring growth rates (cid:38) . t − .The computational time for PEnGUIn is about 12 minutes perorbit on a single GTX-Titan graphics card.
Semi-analytic Method
Equation 16 constitutes an eigenvalue problem, where η m is the eigenfunction and ω is the eigenvalue. To solve thisproblem, we develop a code that directly integrates the di ff er-ential equations, iterates for the correct boundary conditions,and optimizes to find the eigenvalues. The complexity of thiscode is mainly to overcome the di ffi culty imposed by the in-tegral in Equation 16, which e ff ectively raises the order of thedi ff erential equation. The details are documented in AppendixA.Despite the fact that it solves the linearized equations, oursemi-analytic method in fact requires a much longer compu-tational time than simulations using PEnGUIn . Due to limitedresources, initially we only apply it to five sets of parameters.Table 1 contains a list of these sets. One major advantage ofthis method is that it does not have a limit to how slow of agrowth rate can be detected, so we also apply it to all caseswhere simulations do not detect any modal growth. Amongthem, we find a positive growth rate for one case. It is also This case has ( β, c s ) = (0 , . β =
0, the modal growth is purely listed in Table 1, making a total of six sets of parameters. RESULTS
The sets of parameters we consider are β = { , . , . , . , . , . , . } and c s = { . , . , . , . , . } . All 30 combinations of thesevalues are simulated, but only a select few are solved with oursemi-analytic method (see Table 1). The growth rates foundby our two independent approaches agree to ∼ Linear Modes
We find clear growth of asymmetric modes for all caseswith β larger than a certain threshold value that is weakly de-pendent on c s across our parameter space. For most of ourchosen c s values, modal growth is only detected when β ≥ . c s ∼ .
02, where this threshold rises to β ≥ . β and smaller c s , because β measures thestrength of radiation pressure while c s is a source of resis-tance to external forcing. In general, we do find the growthrate to increase with β and decrease with c s , but with obviousexceptions.In Figure 7 we divide our parameter space into three re-gions: regions I and II where modal growth is driven by ra-diation pressure, and region III where it is mainly driven byhydrodynamical e ff ects. In regions I and II, growth rate scalesroughly linearly with β for any given c s , a trend that can bemore easily seen in Figure 8. This is consistent with our ex-pectations.Figure 9, however, reveals a more complicated aspect ofIRI. Disregarding the β = . ≤ c s ≤ .
06 for any given β . Once c s goes below 0 . ff erent trends depending on the value of β : thegrowth rate increases as c s decreases for β ≥ .
2, but for asmaller β the trend flattens or even begins to drop. This com-plex behavior may relate to how sound waves and IRI modescouple. While sound waves have a length scale h , IRI modesare mainly restricted by the sharpness of the transition region,which has a length scale ∆ r . Our results suggest that the cou-pling is weak when h ∼ ∆ r , and becomes much stronger as h decreases, allowing the transition region to accommodate thefull wavelength of the longest wave.Region III is where radiation pressure becomes a smallere ff ect than gas pressure. For c s ≥ .
05, or equivalently, h ( r ) ≥ ∆ r , we detect modal growth even in the absence ofany radiation pressure. In fact, when β =
0, our disk modelis similar to the ”homentropic step jump” model used by Liet al. (2000) to study RWI (see their Figure 2). Their Figure11 shows that for a pressure jump with a width ∆ r = . c s (cid:38) . . Our results are con- hydrodynamical and unrelated to IRI. See Section 5.1 for further discussions. There are a few small di ff erences between the disk model used by Li et al. Fung & Artymowicz F ig . 5.— Fastest growing modes extracted from simulations through Fourier decomposition. Color shows the surface density normalized to the peak of eachmode. On the left is an m =
18 mode from ( β, c s ) = (0 . , . m = β, c s ) = (0 . , . m = β, c s ) = (0 , . ig . 6.— The fastest growing modes directly computed using our semi-analytic method for the same parameters listed in Figure 5. sistent with their findings.The division between IRI and RWI is clear to us because thetwo mechanisms appear to destructively interfere with eachother. For c s = .
06, there is a clear drop in growth rate from β = β = .
05 before it rises again (see Figure 8). Sim-ilarly for c s = .
05, we do not detect any modal growth at β = .
05 even though it is detected at both β = β = . c s and β to haveopposing e ff ects on the epicyclic frequency κ . In Figure 10we see that gas pressure lowers κ near r = r , while β raisesit. Two e ff ects roughly cancel when β = .
05. It is unclearwhether this is coincidental or not.This dividing line may not remain at β = .
05 for a di ff er-ent 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 isunclear to us whether this dividing line will move to a higheror lower β . Additionally, a high c s can push κ below zeroand trigger Rayleigh instability which further complicates thematter. For our disk model this limit is at c s ∼ .
07. Since thefocus of this paper is to characterize IRI, we defer the thor-ough investigation on the interactions between IRI and otherforms of instability to a future study.Other than the growth rate, we also find other general trendsabout the linear modes. For an increasing β or decreasing c s , (2000) and ours. For example, they use an adiabatic equation of state withan adiabatic index of 5 /
3, and their pressure jump is modeled with a di ff erentformula (compare their Equation 3 to our Equation 26). We consider thesedi ff erences insignificant. the azimuthal mode number m of the fastest growing mode in-creases. The dependence on β is particularly pronounced. Inthe most extreme case, the fastest growing mode is m = β, c s ) = (0 . , . β = . m = ∼ m =
18. See Table 1 for more examples. Similarly, we findthat the radial extent of the fastest growing mode becomesmore confined as β increases and c s decreases, as can be seenin Figures 5 and 6. It is therefore empirically apparent that ahigher β encourages the growth of a shorter wavelength mode.While our theory does not make any predictions about whichmode grows the fastest, in hindsight this result is not surpris-ing because the radial motion of an IRI mode must be drivenby radiation pressure, so a stronger perturbing force shouldgenerate a faster radial motion, and therefore a higher fre-quency wave.Another trend is that for an increasing β or decreasing c s ,the corotation radius decreases. This is in accordance withour prediction in Section 2.2. The dependence is weak butnoticeable (see Table 1). Curiously, Figure 5 and 6 show thatthe peak location of each mode is relatively insensitive to vari-ations in both β and c s . Consequently, these peaks generallydo not coincide with their corotation radii.For the bulk of this work we do not explicitly vary ∆ r as afree parameter. Since our choice of ∆ r = .
05 is arbitrary, itis useful to find out to what extent our results would changefor a di ff erent ∆ r . Setting c s = .
04 and ∆ r = .
1, we find therradiation Instability at the Inner Edges of Accretion Disks 7 F ig . 7.— Growth rate of the fastest growing mode as a function of β and c s .The black region is where a positive growth rate is not found with both of ourapproaches. Regions I and II are where IRI operates, while region III seesthe purely hydrodynamical RWI. In the nonlinear phase, clumping occurs inregion I, where local surface density is enhanced by at least a factor of two,often much higher. β I m ( ω ) [ t dyn - ] c s = 0.06c s = 0.05c s = 0.04c s = 0.03c s = 0.02 F ig . 8.— Growth rate of the fastest growing mode as a function of β . The( β, c s ) = (0 , .
05) point is disconnected because no modal growth is detectedat ( β, c s ) = (0 . , . threshold for modal growth becomes β ≥ .
25, roughly twiceas large as when ∆ r = .
05. This suggests that the thresholdvalue scales roughly linearly with ∆ r for the parameter spacewe considered. Nonlinear Evolution
Nonlinear evolution is what separates region I from regionII in Figure 7. Figures 11 and 12 shows the simulation snap-shots for the same two sets of parameters in the left and mid-dle panel of Figure 5 and 6. The left panel, which belongsto region I with ( β, c s ) = (0 . , . Region III is omitted from the discussion to maintain focus on IRI. Werefer the reader to the literature, e.g. Li et al. (2001); Meheut et al. (2012);Lin (2013), for detailed studies on the nonlinear evolution of RWI. s I m ( ω ) [ t dyn - ] β = 0.30 β = 0.25 β = 0.20 β = 0.15 β = 0.10 β = 0.05 β = 0 F ig . 9.— Growth rate of the fastest growing mode as a function of c s . κ / Ω k c s = 0.06 ; β = 0c s = 0.06 ; β = 0.05c s = 0.06 ; β = 0.10 F ig . 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 = very high surface density, exceeding 10 times Σ d , the initialsurface density of the disk defined in Equation 26. This typeof ”clumping”, which we define as a detection of Σ > Σ d anywhere in the disk, is characteristic of region I. Note thatclumping is not a necessary product of IRI, since region II isalso driven by IRI. The transition from region II to I is rapid,in the sense that as we move toward the upper left corner ofFigure 7, the highest local surface density quickly rises to afew tens of Σ d .To compare the two regions in detail, we use the right panelof Figure 11 as a typical case for region II. It shows a signif-icant widening of the edge by a factor of ∼
3. Vortices areformed along the edge and they create mild local enhance-ments in density. Their structure is complex as they typicallylaunch two sets of spiral arms instead of one. The numberof vortices is initially equal to the mode number of the fastestgrowing linear mode, but as they interact with each other, theycan occasionally merge. It is unclear how many will remainin the long run since our simulations only last for 100 orbitsat most.In comparison to region II, the clumping in region I createsa much di ff erent, almost violent, nonlinear evolution. Theclumps are very sharp features, with jumps over three ordersof magnitude in density while their sizes are merely ∼ . r .They are constantly formed and destroyed by disk shear overa dynamical timescale. The destroyed clumps form high den- Fung & Artymowiczsity streams that are also visible in the left panel of Figure11. Another consequence of this clumping is that by concen-trating a large amount of matter in a small region, radiationis able to penetrate further into the disk and push the edge ofthe disk to a higher orbit. In the same figure one can see thatthe edge of the disk is shifted to ∼ . r . We speculate thatthe di ff erence between regions I and II is due to the influenceof gas pressure. Higher gas pressure results in vortex forma-tion more similar to the purely hydrodynamical RWI, and asgas pressure becomes weaker compared to radiation, sharperfeatures are created. ˜ τ m and the Instability Criterion Revisited In our derivation for the instability criterion in Section 2.1,we propose the crude assumption ˜ τ m ∼ τ . Using our semi-analytic method to obtain solutions for η m , we are able toevaluate ˜ τ m explicitly. For all IRI modes we have solvedsemi-analytically, we find ˜ τ m /τ > r = { r − ∆ r , r } , but it is never an order of magnitude above unity.For example, Figure 13 plots ˜ τ m /τ for the m =
18 mode with( β, c s ) = (0 . , . . ≤ r ≤ .
02, Thereal part of ˜ τ m /τ is within 1 to 6, while the imaginary part isclose to vanishing. Using this empirical result we can rewriteEquation 22 as: q β ≈ β e − τ (cid:32) d ln (cid:2) r β (cid:3) d ln r + f τ d ln [ r R ]d ln r (cid:33) , (29)where we substitute ˜ τ m for f τ , and f > f =
3, Figure 14 plots q β for a fewdi ff erent β . This choice of f puts the threshold for instability( q β >
1) at around β = .
1, similar to our empirical results. q β becomes negative as soon as τ > r > r = R quicklydrives its gradient negative. This implies that the instabilitymust originate from the τ < r , asshown in Table 1. CONCLUSIONS AND DISCUSSIONS
We demonstrated that IRI can operate at an inner disk edgewhere there is a transition from being radially transparent toopaque. A local criterion for axisymmetric instability was de-rived (Equation 22). For our given disk model we computedthe 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 from10 − to 10 t − (Figure 7). The fastest rates were found forthe largest β and smallest c s . We empirically determined thatthe threshold for IRI is β ∼ . ∆ r = .
05, with a weakdependence on c s . For a wider edge, ∆ r = .
1, this thresh-old rises to β ∼ .
25. We note that this implies the thresholdcan be lowered by reducing ∆ r ; however, at the same time c s must also be lowered for IRI to dominate over other forms ofinstability that may be triggered by the sharpness of the edge,such as RWI and Rayleigh instability. We employed two in-dependent approaches to obtain the growth rates of the linearmodes: simulating the disks numerically using PEnGUIn , andsolving the linearized equations semi-analytically. Their ex-cellent 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 timesthe local surface density enhancements in the nonlinear evo-lution of IRI.
Connection to Physical Disks
Our disk model is inspired by transitional disks (e.g. Cal-vet et al. 2005; Espaillat et al. 2007, 2008; Andrews et al.2011). The inner edges of these disks are currently unre-solved by observation, but theoretical work has shown thatthe sharpness of disk edges created by X-ray photoevapora-tion (e.g. Owen et al. 2010) is similar to that described by ourEquation 26 with ∆ r = .
05 (compare our Figure 2 to Fig-ure 2 of Owen et al. (2013)). If a transitional disk undergoesIRI, the asymmetric structure at the inner edge will create anazimuthal variation in shadowing. Flaherty and Muzerolle(2010) showed that this can lead to a significant variation indisk emission. Indeed, some variability in the infrared emis-sion 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 circumstel-lar disks. AGN accretion disks, for example, can be subjectedto IRI if there are any sharp jumps in density and / or opacity,such as the inner edges of the board-line regions. IRI can po-tentially generate the stochastic asymmetry, which is used toexplain the variability in the double-peaked Balmer emissionlines in radio-loud AGNs (Flohic and Eracleous 2008). Wenote that the dynamics in AGN accretion disks are consid-erably more complicated since they do not have a point-likelight source. Implications of ”Clumping”
The ”clumping” found in a part of our parameter space (Fig-ure 7) opens new possibilities for IRI. For instance, very highdensity regions in protoplanetary disks may be favorable envi-ronments for the formation of planetary cores. The density ofindividual clumps may even become high enough to triggergravitational instability at the inner edges of massive disks.One should be cautious to interpret the enhancement factorsreported as realistic, however, since it is only one disk modelthat we have studied.The clumping also leads to a possibility of preventing in-ward dust migration. Dominik and Dullemond (2011) demon-strated that while radiation pressure can initially push dustoutward and form a dust wall, the wall eventually succumbsto the global accretion flow and migrates inward. If this wallbecomes unstable due to IRI, clumping can occur, e ff ectivelycreating ”leakage” within the wall, allowing radiation to pushdust further back. The true behavior of these dust walls is im-portant to understand disks where inner clearings have beenobserved, such as transitional disks. Dynamical interactionsbetween radiation, dust, and gas must be considered for thiskind of study. Outlook
There are three main aspects of our model that we feelwould benefit greatly from a more realistic treatment. First,our model ignores the vertical dimension. A notable di ff er-ence from 2D to 3D is that the location of the inner edge ofa disk, defined as the τ = ∼ h . One possibleconsequence is that IRI would generate a vertical circulationat the inner edge, which would dilute the opacity in the mid-plane and allow radiation pressure to penetrate further into thedisk. Additionally, in a flared disk, radiation pressure is ex-erted on the photosphere of the entire disk rather than just theinner edge. On the other hand, because of dust settling, weexpect the value of β in the photosphere to be smaller than therradiation Instability at the Inner Edges of Accretion Disks 9 F ig . 11.— Snapshots of our simulations for ( β, c s ) = (0 . , .
02) on the left and ( β, c s ) = (0 . , .
05) on the right, taken at t =
100 orbits. Surface density is shownin logarithmic scale. The simulation on the left, belonging to region I of Figure 7, shows very high local surface density, an e ff ect we describe as ”clumping”.On the right, belonging to region II of Figure 7, shows 6 vortices with di ff erent orbital frequencies but all lining up near r = . ∼ .
2. Each of these vorticeslaunches two pairs of spiral arms.F ig . 12.— Cartesian view of Figure 11. midplane, making it even more di ffi cult to reach the β ≈ . τ = ff ect of radiation pressure,and in the optically thick disk, dust migrates inward due togas drag. This behavior of dust is described in Section 3 ofTakeuchi and 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 heateddirectly by irradiation, to passively by the irradiated atmo-sphere. Consequently the midplane temperature should bedecreasing across the disk edge. This is not captured by ourglobally isothermal assumption. Additionally, the clumps wefind in some of our nonlinear results are su ffi ciently dense thatthey are optically thick. With our isothermal treatment, theyremain the same temperature as their surroundings, while in0 Fung & Artymowicz r -10-8-6-4-20246810 τ m / τ realimaginary ~ F ig . 13.— ˜ τ m /τ for the m =
18 mode with ( β, c s ) = (0 . , . r ≈ { . , . } , the approximation ˜ τ m ∼ τ is accurate towithin order unity. r q β β = 0.05 β = 0.10 β = 0.15 β = 0.20 β = 0.25 β = 0.30 F ig . 14.— q β from Equation 29 for di ff erent values of β . We choose f = truth these clumps should be capable of shielding themselvesfrom irradiation and creating a non-trivial internal tempera-ture structure. Whether this is an e ff ect that aids or inhibitstheir formation and survival requires future investigation. ACKNOWLEDGMENTS
We thank Yanqin Wu, Chris Matzner, and Eugene Chi-ang for helpful feedback that substantially improved thismanuscript. We also thank James Owen for insightful discus-sions. JF owes his gratitude to the Queen Elizabeth II Grad-uate Scholarship in Science and Technology. We gratefullyacknowledge support from the Discovery Grant by the Natu-ral Sciences and Engineering Research Council of Canada.
APPENDIXNUMERICAL 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 benumerically integrated in the direction of increasing r because of the integral in the fourth term. In principle, it is possible tosimply 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 integrateEquation 16 simultaneously from the inner boundary outward, and the outer boundary inward, and find the ω that results in amatch of the two functions at some intermediate radius r mid . To accomplish this, we first define y m ≡ (cid:90) r η m c d τ dr (cid:48) dr (cid:48) , (A1)and then di ff erentiate Equation 16 with respect to r : ∂ y m ∂ r + a (cid:48) ( r ) ∂ y m ∂ r + b (cid:48) ( r ) ∂ y m ∂ r + c (cid:48) ( r ) y m = , (A2)where a (cid:48) ≡ a − d ln gdr , b (cid:48) ≡ b − a d ln gdr + (cid:32) d ln gdr (cid:33) − g d gdr , c (cid:48) ≡ cg , g ≡ c d τ dr . Thus we can now numerically integrate Equation A2 in both directions, and recover η m from y m . The boundary conditions canrradiation Instability at the Inner Edges of Accretion Disks 11be approximated using the WKB method. The WKB form for y m is y m = R ( r ) e i (cid:82) r kdr (cid:48) , (A3) ∂ y m ∂ r (cid:39) iky m , (A4)where R ( r ) is a slowly varying function and k is the complex wave number that satisfies | kr | (cid:29)
1. Substituting Equation A3 andEquation A4 into Equation A2 we get the following algebraic equation for k : k − ia (cid:48) k − b (cid:48) k + ic (cid:48) = . (A5)The three solutions of Equation A5 correspond to the inward traveling (Re( k ) < k ) > | kr | (cid:28)
1, e ff ectively rendering y m a constant,which violates the approximation of a tightly winding wave. To accommodate for this solution, we generalize Equation A3 toallow for a constant o ff set: y m = R ( r ) e i (cid:82) r kdr (cid:48) + C . (A6)Substituting this into Equation A2, we obtain: k − ia (cid:48) k − b (cid:48) k + ic (cid:48) (cid:32) y m y m − C (cid:33) = . (A7)In the optically thin and thick limits, c (cid:48) becomes arbitrarily small, and since the | kr | (cid:28) ff set C , the last term can be dropped, giving back the usual quadratic form: k − ia (cid:48) k − b (cid:48) = , (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 . Forclarity, we will denote variables associated with the solution integrated from the inner boundary with the subscript ”in”, and thethose 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 inwardof 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 thefollowing iterative formulas: R i + = R i out d y im , in dr d y im , out dr − , (A9) C i + = C i out + y im , in − y im , out , (A10)where i is the current iterative step, the y m and its derivatives are evaluated at r mid . Convergence typically requires tens or evenhundreds of iterations, which is the primary reason for the large amount of computational time required for this method. Lastly,we find the eigenvalue ω by minimizing the following function, evaluated at r mid : f = Re (cid:16) dy m , out dr (cid:17) − Re (cid:16) dy m , in dr (cid:17) max (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) Re (cid:16) dy m , out dr (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12)(cid:12) Re (cid:16) dy m , in dr (cid:17)(cid:12)(cid:12)(cid:12)(cid:12)(cid:21) + Im (cid:16) dy m , out dr (cid:17) − Im (cid:16) dy m , in dr (cid:17) max (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) Im (cid:16) dy m , out dr (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12)(cid:12) Im (cid:16) dy m , in dr (cid:17)(cid:12)(cid:12)(cid:12)(cid:12)(cid:21) . (A11)We use an eighth-order Runge-Kutta method with adaptive step-size control for the numerical integration. We set r mid = r in = .
3, and the outer at r out =
4. Minimizing f is also very time consuming because we employ arandom sampling method: first we bracket the minimum within a range of likely values for the real and imaginary part of ω , thenwe 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 thegradient of f , because of the numerous local minima that exist. REFERENCESS. A. Balbus and J. F. Hawley, Reviews of Modern Physics , 1 (1998).D. N. C. Lin and J. E. Pringle, MNRAS , 607 (1987).C. F. Gammie, ApJ , 174 (2001), arXiv:astro-ph / , 721 (1984).J. C. B. Papaloizou and J. E. Pringle, MNRAS , 799 (1985).J. C. B. Papaloizou and J. E. Pringle, MNRAS , 267 (1987).P. Goldreich, J. Goodman, and R. Narayan, MNRAS , 339 (1986).R. V. E. Lovelace, H. Li, S. A. Colgate, and A. F. Nelson, ApJ , 805(1999), arXiv:astro-ph / , 399 (1998). N. Higginbottom, D. Proga, C. Knigge, K. S. Long, J. H. Matthews, andS. A. Sim, ArXiv e-prints (2014), 1402.1849.P. Thebault, Q. Kral, and J.-C. Augereau, A&A , A16 (2014), 1310.1584.E. Chiang and R. Murray-Clay, Nature Physics , 604 (2007), 0706.1241.C. Dominik and C. P. Dullemond, A&A , A101 (2011), 1105.4459.F. Kirchschlager and S. Wolf, A&A , A54 (2013), 1302.5275.P. Goldreich and S. Tremaine, ApJ , 857 (1979).J. Fung, J.-M. Shi, and E. Chiang, ApJ , 88 (2014), 1310.0156.P. Colella and P. R. Woodward, Journal of Computational Physics , 174(1984).J. M. Blondin and E. A. Lufkin, ApJS , 589 (1993). H. Li, J. M. Finn, R. V. E. Lovelace, and S. A. Colgate, ApJ , 1023(2000), arXiv:astro-ph / ff , and R. Liska, ApJ , 874 (2001),astro-ph / , A134 (2012),1208.4947.M.-K. Lin, ApJ , 84 (2013), 1209.0470.N. Calvet, P. D’Alessio, D. M. Watson, R. Franco-Hern´andez, E. Furlan,J. Green, P. M. Sutter, W. J. Forrest, L. Hartmann, K. I. Uchida, et al.,ApJL , L185 (2005).C. Espaillat, N. Calvet, P. D’Alessio, E. Bergin, L. Hartmann, D. Watson,E. Furlan, J. Najita, W. Forrest, M. McClure, et al., ApJL , L111(2007), 0707.0019.C. Espaillat, J. Muzerolle, J. Hern´andez, C. Brice˜no, N. Calvet,P. D’Alessio, M. McClure, D. M. Watson, L. Hartmann, and B. Sargent,ApJL , L145 (2008), 0810.4575. S. M. Andrews, D. J. Wilner, C. Espaillat, A. M. Hughes, C. P. Dullemond,M. K. McClure, C. Qi, and J. M. Brown, ApJ , 42 (2011), 1103.0284.J. E. Owen, B. Ercolano, C. J. Clarke, and R. D. Alexander, MNRAS ,1415 (2010), 0909.4309.J. E. Owen, M. Hudoba de Badyn, C. J. Clarke, and L. Robins, MNRAS , 1430 (2013), 1309.0508.K. M. Flaherty and J. Muzerolle, ApJ , 1733 (2010), 1007.1249.J. Muzerolle, K. Flaherty, Z. Balog, E. Furlan, P. S. Smith, L. Allen,N. Calvet, P. D’Alessio, S. T. Megeath, A. Muench, et al., ApJL , L15(2009), 0909.5201.K. M. Flaherty, J. Muzerolle, G. Rieke, R. Gutermuth, Z. Balog, W. Herbst,S. T. Megeath, and M. Kun, ApJ , 83 (2011), 1103.0781.C. Espaillat, E. Furlan, P. D’Alessio, B. Sargent, E. Nagel, N. Calvet, D. M.Watson, and J. Muzerolle, ApJ , 49 (2011), 1012.3500.H. M. L. G. Flohic and M. Eracleous, ApJ , 138 (2008), 0806.0163.T. Takeuchi and P. Artymowicz, ApJ , 990 (2001), astro-ph //