# The Irradiation Instability of Protoplanetary Disks

Yanqin Wu  Department of Astronomy and Astrophysics, University of Toronto,
Toronto, ON M5S 3H4, Canada  Yoram Lithwick  Dept. of Physics and Astronomy,
Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208 & Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA)

###### Abstract

The temperature in most parts of a protoplanetary disk is determined by
irradiation from the central star. Numerical experiments of Watanabe & Lin
(2008) suggested that such disks, also called ‘passive disks’, suffer from a
thermal instability. Here, we use analytical and numerical tools to elucidate
the nature of this instability. We find that it is related to the flaring of
the optical surface, the layer at which starlight is intercepted by the disk.
Whenever a disk annulus is perturbed thermally and acquires a larger scale
height, disk flaring becomes steeper in the inner part, and flatter in the
outer part. Starlight now shines more overhead for the inner part and so can
penetrate into deeper layers; conversely, it is absorbed more shallowly in the
outer part. These geometric changes allow the annulus to intercept more
starlight, and the perturbation grows. We call this the irradiation
instability. It requires only ingredients known to exist in realistic disks,
and operates best in parts that are both optically thick and geometrically
thin (inside $30$AU, but can extend to further reaches when, e.g., dust
settling is considered). An unstable disk develops travelling thermal waves
that reach order-unity in amplitude. In thermal radiation, such a disk should
appear as a series of bright rings interleaved with dark shadowed gaps, while
in scattered light it resembles a moving staircase. Depending on the gas and
dust responses, this instability could lead to a wide range of consequences,
such as ALMA rings and gaps, dust traps, vertical circulation, vortices and
turbulence.

##  1 Introduction

Currently, the main bottleneck for understanding planet formation lies in an
incomplete knowledge of the protoplanetary disk. In this work, we study the
dynamics of passive disks, i.e., disks where stellar irradiation dominates the
energetics. These include almost all parts of the disk except perhaps for the
innermost region. We show that passive disks suffer from an instability, and
this could have a broad range of observational and theoretical ramifications.

###  1.1 Motivation

I. Gaps and Rings: Recent observations made using the Atacama Large Millimetre
Array (ALMA) have shown that typical disks are not the smooth power-laws
beloved by theorists. Rather, bright rings and dark gaps are ubiquitous, on
scales from tens to hundreds of AU (e.g., ALMA Partnership et al., 2015;
Andrews, 2020; Huang et al., 2018). At the moment, these rings and gaps are
most commonly attributed to the effects of unseen planets (e.g., Baruteau et
al., 2014; Dong et al., 2015; Dipierro et al., 2016; Bae et al., 2017; Dong et
al., 2017; Zhang et al., 2018). And in a few cases, there is strong kinematic
evidence for planets, such as in the gap of HD 163296 (Pinte et al., 2018).
But the near-ubiquity of gaps and rings (Long et al., 2018; Nielsen et al.,
2019) is at tension with the paucity of large-mass planets at these distances,
as suggested by direct imaging surveys (Bowler, 2016) and microlensing surveys
(Suzuki et al., 2016; Gaudi, 2021). While the planet hypothesis is difficult
to exclude, given its large number of free parameters (such as planet mass and
orbit, disk viscosity, and history) it is worthwhile to consider whether
planets are the causes for the gaps and rings, or instead the products of such
features.

Many alternative scenarios have been proposed to explain these features,
including dust-drift-driven viscous ring instability (Wünsch et al., 2005;
Dullemond & Penzlin, 2018), secular gravitational instabilities in the dust
(Takahashi & Inutsuka, 2014), dead zones (Flock et al., 2015), snow lines
(Okuzumi et al., 2016), MHD wind-driven structures (Bai, 2014; Riols et al.,
2020), and an eccentric disk instability (Li et al., 2021). Lastly,
Siebenmorgen & Heymann (2012) and Ueda et al. (2019) propose that they may be
triggered by an instability found in Monte Carlo simulations of irradiated
disks. This last proposal may be closely related to the irradiation
instability considered in this paper.

These features motivate us to study the stability of passive disks. An
inherent instability could naturally explain the multiple gaps and rings in a
given disk, without invoking an arbitrary number of planets. It would also
directly impact the formation and migration of planets.

II. Dust Wafting and Migration: Radiation in disks is controlled by dust.
Micron-sized grains absorb and scatter starlight, while larger ($\sim$mm-
sized) grains provide the bulk of the opacity for the thermalized radiation.
Conversely, the distribution of dust grains is strongly influenced by gas
dynamics.

Both the spectral energy distribution (e.g., Chiang et al., 2001; Furlan et
al., 2006; Woitke et al., 2019) and scattered light images (e.g., Avenhaus et
al., 2018) of protoplanetary disks suggest that micron-sized grains must
reside at least a couple gas scale heights above the midplane. Although they
are tightly coupled to the gas, vertical settling high up in the atmosphere is
rapid (Dullemond & Dominik, 2004a), even more so when dust coagulation occurs
(Dullemond & Dominik, 2005). These small grains must be replenished, by
fragmentation and/or up-draft. Both can be provided by gas turbulence or
circulation.

The mm-sized dust, on the other hand, is much more weakly coupled to the gas,
and so settles closer to the midplane. Without any local pressure maxima,
these grains can drift for large distances over the disk’s lifetime
(Weidenschilling, 1977). Pressure maxima, on the other hand, can halt this
migration and trap drifting particles (Whipple, 1972). ALMA images
tantalizingly suggest that the observed rings and gaps are the smoking gun for
particle traps (e.g., Andrews & Birnstiel, 2018; Dullemond et al., 2018).

A disk instability could therefore have important implications for the dust
behaviour in protoplanetary disks. It could generate turbulence or meridional
flows that waft up the micron-sized grains. And it may also produce pressure
maxima, leading to vortices or axisymmetric rings, natural barriers for
inwardly drifting dust and welcoming cradles for planet formation.

III. Turbulence and Accretion: A crucial open question in the study of
protoplanetary disks is why disks accrete (Lynden-Bell & Pringle, 1974). It is
unclear if these disks are turbulent, and whether such turbulence can
transport enough angular momentum to disperse disks in a few million years
(see review by Klahr et al., 2018). Until recently, the most promising
mechanism for generating turbulence was the magnetorotational instability
(Balbus & Hawley, 1998). But such MHD proposals require a sufficient degree of
ionization that the disk gas can couple effectively to the magnetic field.
That is hard to achieve, particularly for disks that are very dusty (e.g.,
Bai, 2015; Simon et al., 2015). A plethora of other instabilities have been
investigated (as reviewed in Klahr et al., 2018). One example is the vertical
shear instability (Urpin & Brandenburg, 1998; Lin & Youdin, 2015), which may
produce sufficiently strong turbulence for accretion ($\alpha\sim 10^{-4}$,
where $\alpha$ is the Shakura-Sunyaev parameter), but at large distances
(Flock et al., 2020). This $\alpha$-value appears consistent with some upper
limits placed at various disk locales (e.g., $\alpha\lesssim 7\times 10^{-3}$,
Flaherty et al., 2018). Alternatively, disks may also disappear via a non-
turbulent mechanism such as disk winds (Bai & Stone, 2013).

As the situation remains murky, we are motivated to look for a robust
instability that could generate turbulence or fluid circulation, and
ultimately drive accretion.

###  1.2 Prior Work on Passive Disks and Their Stability

For a protoplanetary disk accreting at a typical rate (say, $\dot{M}\sim
10^{-8}M_{\odot}$/yr), disk heating is dominated by stellar irradiation except
inward of $\sim 1$AU (see, e.g., D’Alessio et al., 1998). So, most of the disk
is should be regarded as being ‘passively heated.’ Kenyon & Hartmann (1987)
showed that passive disks can account for the far-IR excesses of IRAS disks,
provided they are flared. The predicted flaring morphology was confirmed by
Hubble images (e.g., Burrows et al., 1996).

Chiang & Goldreich (1997, hereafter CG97) set forth a simple model for such
disks. Optical light from the star is absorbed high up in the disk by small
grains. And radiation from these grains illuminates the disk midplane. In
thermal and hydrostatic equilibrium, flared passive disks take a simple
analytic form: $\frac{h}{r}\propto r^{2/7}$, where $h$ is the vertical scale
height and $r$ the cylindrical radius (Kusaka et al., 1970; Cunningham, 1976;
Chiang & Goldreich, 1997). Passive disk models have been very successful in
explaining the spectral energy distributions of protoplanetary disks, provided
one also accounts for some vertical settling of the dust (Chiang et al., 2001;
Dullemond & Dominik, 2004a; D’Alessio et al., 2006).

The stability of passive disks was first investigated by Cunningham (1976) for
disks around black holes. That work was extended to protoplanetary disks by
D’Alessio et al. (1999). Their simple analysis showed that such disks are
stable: thermal perturbations propagate inwards, and damp along the way.

The equilibrium solution for a stable disk should be easily obtained by
iteration. Yet mysteriously, such attempts are often plagued by convergence
issues (e.g., Dullemond & Dominik, 2004b; Min et al., 2009; Siebenmorgen &
Heymann, 2012; Wang & Goodman, 2017; Ueda et al., 2019). Using Monte Carlo
radiative transfer codes to describe the radiation effects more accurately
than CG97, these authors iteratively solved the equations of hydrostatic
equilibrium and thermal equilibrium. They often find no convergence,
particularly for disks with realistically large dust surface densities. With
successive iterations, new waves appear at large radii, and propagate inward
with order-unity amplitudes (see, e.g., Fig. 7 of Ueda et al., 2019). It is
unclear if such behaviour is generic in physical disks, or if it is an
artificial instability introduced by the iteration procedure. In any case,
this issue hampers further study of realistic disks.

The work by Watanabe & Lin (2008), though receiving little attention (see Ueda
et al., 2021, for a modern reincarnation), raises an interesting possibility.
Using 1-D simplified radiative transfer, half-way in complexity between CG97
and a Monte Carlo code, they found that passive disks are unstable. They
argued that the instability is likely related to changes in the optical
surface. Their numerical experiments showed that the unstable disk develops
inward travelling thermal-waves that reach order unity amplitudes.

In this work, we further elucidate the origin of this behaviour, which we call
the ‘irradiation instability.’ Using both analytical and numerical tools
(including the radiative transfer code RADMC-3D), we demonstrate that the
instability is genuine, not numerical – although using RADMC-3D (Dullemond et
al., 2012, henceforth called RADMC) to iterate is risky unless guided by
analytical insights. We derive the conditions for such an instability and
argue that they should be prevalent in observed disks.

Paper Overview: In §2 we present a cartoon view of the instability. In §3–6 we
do the math. Because of the complexities involved with solving the radiative
transfer problem, we develop three models that are successively more complex,
and more realistic. Readers who prefer to skip the technical details may
proceed to §7 where we summarize the main results. We end with an extensive
discussion of the assumptions (§8) and a brief introduction of things to come
(§9).

##  2 The Irradiation Instability – Cartoon Version

###  2.1 Passive Disks in Equilibrium

Figure 1: Cartoon of a passive disk, in equilibrium (left) and perturbed
(right). The color shading illustrates the dust density relative to that at
midplane. Here, $h$ is the vertical gas scale height, while $H$ (red line) is
the ‘optical surface,’ the height at which radial starlight (green lines) is
intercepted by the disk. The angular extent of the white segment, as seen by
the star, indicates the amount of starlight intercepted by the disk in that
radial zone. As this zone experiences a rise in midplane temperature (right
panel), the amount of starlight it receives increases, as evidenced by the
larger angle subtended by the two green rays. This increase occurs for two
reasons: the rise of the optical surface in proportion to the rise in scale
height (dashed curve), and changes in the penetration depths of the star rays.

We follow the concepts and notations of CG97 to study passive disks. Fig. 1
(left panel) shows a cartoon of a passive disk in equilibrium, with the
background shading representing the density of dust grains, relative to their
midplane density. The star’s optical light is absorbed by dust grains at
altitude $H$, the ‘optical surface,’ that lies a few scale heights above the
midplane. This layer re-radiates half of the luminosity it receives down into
the disk, and the latter then re-emits this energy at longer wavelengths
(thermalized radiation).

The amount of heating a disk receives is determined by the flaring of its
optical surface. In Fig. 1, the relevant star-rays received by a radial zone
of concern are those bound by the two green arrows. When balancing the heating
against blackbody cooling, and insisting on vertical hydrostatic equilibrium,
a flared solution is found for these disks (Kusaka et al., 1970; Chiang &
Goldreich, 1997; Dullemond, 2000).

###  2.2 Passive Disk Perturbed

Fig. 1 illustrates what happens to a disk when a localized thermal
perturbation increases the scale-height $h(r)$. There are two effects. First,
the optical surface rises in proportion to the scale height (dashed curve).
This is what is considered in D’Alessio et al. (1999). The disk intercepts a
bit more stellar flux, but not enough to overcome the extra blackbody cooling
from a now hotter disk. As a result, the perturbation damps away.

But there is another effect. Consider first the inner half of the scale height
perturbation. The slope of the optical surface is increased there, so the
star’s rays shine more directly overhead (i.e., closer to the surface normal),
and can penetrate more deeply into the disk. This effect is analogous to
stellar limb-darkening, but now for absorption. Conversely, at the outer part
of the perturbation, starlight penetrates more shallowly due to the more
tangent slant. As a result, the opening angle between the two green rays is
increased. This means excess heating, and given the right perturbation, it can
overcome the excess cooling and drive an instability.

Interestingly, an extreme manifestation of stronger flaring leading to
enhanced heating is the inner rim of a protoplanetary disk (Dullemond et al.,
2001). In this region, the star’s rays enter the disk almost head-on. The
strong heating puffs up the inner rim into a wall.

In the following sections, we make this picture quantitative.

##  3 Setting the Stage

To study the stability of a passively irradiated disk, we make a number of
simplifying assumptions, both to highlight what we believe are the most
relevant effects, and to facilitate a short treatment. These include:

– opacity is provided only by dust;

– dust traces gas with a constant ratio;

– vertical hydrostatic equilibrium;

– the radial profile of surface density does not vary;

– axisymmetry;

– gas temperature tracks dust temperature;

– star is point-like so all star rays are radial;

– no inner hot rim, which would otherwise cast a shadow.

In addition to these, we also introduce a string of other simplifications,
which we will describe as we go along. In §8, we assess how some of these
assumptions and how varying them may qualitatively affect our results.

###  3.1 The Thermal Equation

We adopt cylindrical co-ordinates ($r$ and $z$), and assume that the disk is
sufficiently thin that we may work to leading order in $z/r$ (the small-angle
approximation). The midplane temperature $T$ is governed by the thermal
equation (obtained in Appendix A under a number of assumptions)

$\frac{3}{8}\,\tau_{\rm mm}c_{p}\Sigma_{\rm
gas}\frac{{\partial}T}{{\partial}t}={\frac{1}{2}F_{\rm irr}[T]}-\sigma_{\rm
SB}T^{4}\ .$ (1)

This describes the rate of thermal energy increase per unit disk area as a
result of the imbalance between starlight heating and black-body cooling.
Here, $c_{p}$ is the specific heat per unit mass, $\Sigma_{\rm gas}$ is the
gas surface density, and $\tau_{\rm mm}$ is the vertical optical depth for
thermal radiation (more details below), which is assumed to be larger than
unity. On the right-hand side, $F_{\rm irr}$ is the stellar flux incident on
the disk’s optical surface. The square brackets denote that $F_{\rm irr}$
depends on the local profile (not just value) of temperature. And the factor
of $1/2$ multiplying $F_{\rm irr}$ arises because half of the incident flux is
re-radiated by grains in this surface to space, without heating the disk
interior.

In writing eq. (1), we have simplified the physics of vertical transport under
the assumption that the timescale of variation is much longer than the thermal
time (see Appendix A). A more accurate treatment should allow for vertical
thermal waves. We argue in Appendix A.3 that this likely does not impact our
results.

As we show in Appendix A, the $T$ that appears throughout equation (1) should
be the midplane temperature, not the surface temperature, provided the factor
of $\tau_{\rm mm}$ is included on the left-hand side.111 The factor of
$\tau_{\rm mm}$ has been incorrectly neglected in the literature (D’Alessio et
al., 1999; Dullemond, 2000; Watanabe & Lin, 2008). But its neglect leads only
to a change in timescale.

The heating flux is (Safronov, 1962; Kusaka et al., 1970; Chiang & Goldreich,
1997)

$\displaystyle F_{\rm irr}=\frac{L_{*}}{4\pi r^{2}}\,\sin\alpha\,,$

where $\alpha$ is the grazing angle of the star-rays relative to the optical
surface. We denote the height of this surface above the midplane by $H$. In
the limit that the disk is thin ($H/r\ll 1$) and the star is point-like, we
can simplify the expression for $\alpha$ to arrive at (CG97, Chiang et al.,
2001),

$\displaystyle F_{\rm irr}{\approx}\frac{L_{*}}{4\pi
r^{2}}\,\left(\frac{d}{d\ln r}\,\frac{H}{r}\right)\ .$ (2)

If $H/r$ decreases with radius, the disk falls into shadow cast by interior
annuli, and one should instead set $F_{\rm irr}\rightarrow 0$. Hidden in this
equation is the implicit assumption that the stellar heating is only
determined by the local gradient of the optical surface. This is reasonable as
long as the variation length scale is much greater than $H$, as is certainly
satisfied by the equilibrium disk. However, this may not be true for a
perturbed disk, an issue we call ‘horizontal averaging’ and return to later.

In hydrostatic equilibrium, the midplane temperature determines the gas scale
height as

$\displaystyle\frac{h}{r}=\frac{c_{s}}{\Omega r}=\left(\frac{k_{B}}{\mu
m_{p}GM_{*}}\right)^{1/2}\sqrt{rT}\ ,$ (3)

where $M_{*}$ is the stellar mass, $m_{p}$ the proton mass, and $\mu$ the mean
molecular weight (henceforth set to $\mu=2.3$). The scale height in turn sets
the dust density field, which we take to be

$\rho_{\rm dust}(r,z)=\frac{\Sigma_{\rm
dust}(r)}{\sqrt{2\pi}h(r)}\,\exp\left[-\frac{z^{2}}{2h(r)^{2}}\right]\,,$ (4)

where $\Sigma_{\rm dust}$ is the vertical column density of dust (both sides
of equator). The dust density controls how far the stellar flux penetrates
into the disk.

In this study, we assume that vertical hydrostatic equilibrium is maintained
instantaneously, or eqs. (3)-(4) remain valid as the disk heats and cools.
This is likely valid in the region where the thermal time is longer than the
orbital time, or inwards of $\sim 50$AU for our fiducial disk (see below).

###  3.2 The Optical Surface

The key ingredient in this problem is the geometry of the optical surface. For
our analytical study (not the RADMC simulations), we consider two approximate
forms for $H$:

  1. 1.

Simplistic Surface: following CG97, we introduce an important quantity $\chi$,
the ratio between the optical surface height and the local gas scale height,

$\displaystyle\chi\equiv{H\over h}\ .$ (5)

The value of $\chi$, for passive disks in equilibrium, depends on dust density
logarithmically (see more below). CG97 found that the value of $\chi$ ranged
from $5$ to $4$ between 3AU and 100AU, for a MMSN-type disk. So at first sight
it seems reasonable to assume that $\chi$ is a constant everywhere, and does
not vary when the disk is perturbed. We call this a ‘simplistic surface.’ It
is adopted by D’Alessio et al. (1999) for their stability analysis.

  2. 2.

Realistic Surface: alternatively, one could self-consistently determine the
optical surface using the definition that the optical depth to the star is
unity. This brings in some algebraic difficulties but, as we show below, is
essential for describing perturbed disks. The name ‘realistic’ is a euphemism
– this approach remains an approximation to reality, which can be addressed
only using radiative transfer codes.

In §4-5, we study passive disks under each of these approximations. The key
result is that under the first approximation the disk is stable, whereas under
the second more accurate one it is unstable. We then perform RADMC simulations
(§6) to demonstrate that true disks also exhibit instability.

###  3.3 Fiducial Disk

We assume the star has solar mass and luminosity, $M_{*}=M_{\odot}$,
$L_{*}=L_{\odot}$. For our fiducial disk, we choose a dust surface density

$\displaystyle\Sigma_{\rm dust}$ $\displaystyle=$ $\displaystyle 20\,r_{\rm
AU}^{-1}\ \ \textrm{g}/{\rm cm}^{2}\,$ (6)

where $r_{\rm AU}=(r/{\rm 1AU})$. This density is similar to the minimum mass
solar nebula at 1AU (Hayashi, 1981), but falls off more gradually with radius,
and so is more consistent with many observed systems (e.g., Cleeves et al.,
2016; van Boekel et al., 2017).

We adopt the dust opacity from Fig. 2 of Woitke et al. (2016), calculated for
a power-law dust mixture with sizes from $0.05\mu m$ to $3$mm, which we
approximate as

$\displaystyle\kappa_{\lambda}$ $\displaystyle=$ $\displaystyle
10^{3}\left({\lambda\over 0.5\mu m}\right)^{-1/2}{\rm
cm}^{2}/\textrm{g(dust)}\ ,$ (7)

for total extinction. In the optical, this opacity falls below that adopted by
CG97 by a factor of 40, as they assume all grains are small ($0.1\mu$m).

In our analytical study, the radiation field is described by fluxes at only
two frequencies: that of starlight and of disk thermal radiation. So the disk
radiative properties can be encapsulated by two vertical optical depths. One
is $\tau_{V}$, the optical depth in the visual band,

$\displaystyle\tau_{\rm V}\equiv{1\over 2}\,{{\kappa_{V}\Sigma_{\rm
dust}}}=10^{4}\,r_{\rm AU}^{-1}\ ,$ (8)

where the factor of $1/2$ indicates integration from the midplane upward. The
other is the optical depth for dust thermal radiation. We name it $\tau_{\rm
mm}$, with the actual wavelength determined by the peak of the local
blackbody. Fig. 2 (bottom panel) shows these two optical depths. Our fiducial
disk remains optically thick to thermal radiation out to $\sim 100$AU.

Figure 2:  Our fiducial disk. The top panel compares the thermal time (red
curve) and the dynamical time (orbital period, green curve); inward of $\sim
50$AU, it may be reasonable to assume that vertical hydrostatic equilibrium is
reached quickly. The bottom panel displays the vertical optical depths for
visual (blue) and for thermal radiation (red). The disk is optically thick
throughout. Here, the disk temperature profile is that of the equilibrium
disk, and $\chi={H\over h}$, both determined in §5.1.

Equation (1) induces us to define an evolutionary timescale, the thermal time,

$\displaystyle t_{\rm th}={3\over 8}{c_{p}\,\Sigma_{\rm gas}\,\tau_{\rm
mm}\over\sigma_{\rm SB}T^{3}}\,.$ (9)

We evaluate this timescale using the equilibrium disk temperature, and a gas-
to-dust ratio $\Sigma_{\rm gas}/\Sigma_{\rm dust}=100$, and plot it in Fig. 2.
It decays outward as $t_{\rm th}\propto r^{-0.93}$, and intersects the disk
dynamical timescale at around $50$AU. Outside this radius, the assumption of
vertical hydrostatic equilibrium fails.

##  4 Stability for Simplistic Surface

We solve the thermal equation under the assumption that $\chi=H/h$ is a
constant in both space and time. While allowing for a spatially varying $\chi$
does not much alter our conclusions, a time-varying $\chi$ will, as we discuss
later. We will present the equilibrium and a simple stability analysis,
neither of which is new to this work. We do so both to connect to the previous
(mis)-understanding that passive disks are stable, and to highlight the
differences between such a treatment and a more correct one (§5).

###  4.1 Equilibrium Disk

The equilibrium solution has been presented by Kusaka et al. (1970), CG97, and
Dullemond (2000). The thermal equation reads

$\displaystyle{L_{*}\over 8\pi\sigma_{\rm SB}}{d\over d\ln r}{H\over
r}=r^{2}T^{4}\ .$ (10)

It may be integrated after setting $H=\chi h$ . With $\chi$ being a constant,
and $h$ determined by hydrostatic equilibrium (eq. 3), this equation has a
power-law solution,222 The general solution for $h/r$ is given by eq. 11
multiplied by the factor $\left(1+ar^{2}\right)^{-1/7}$, where $a$ is an
integration constant (Dullemond, 2000). Here we focus only on the power-law
part, as it is relevant for most of the disk.

$\displaystyle{h\over r}=0.02\,\,\chi^{1/7}\,r_{\rm
AU}^{2/7}\times\left({L_{*}\over
L_{\odot}}\right)^{1/7}\times\left({M_{\odot}\over M_{*}}\right)^{4/7}\ ,$
(11)

and

$\displaystyle T=90{\rm K}\,\chi^{2/7}\,r_{\rm
AU}^{-3/7}\times\left({L_{*}\over
L_{\odot}}\right)^{2/7}\times\left({M_{\odot}\over M_{*}}\right)^{1/7}\ .$
(12)

These power-laws reproduce eqs. (14) in CG97, albeit with slightly different
normalizations.

As described below (§5.1), $\chi$ can be solved for self-consistently. It
varies slowly for a disk in equilibrium, so eqs. (11)–(12) remain
approximately correct. In particular, for our fiducial disk, $\chi$ may be
adequately fitted as $\chi\approx 4.5r_{\rm AU}^{-.075}$, and we estimate the
optical surface to lie at

${H\over r}=\,\chi{h\over r}\,\approx 0.11\,r_{\rm AU}^{0.21}\,.$ (13)

###  4.2 Linear Perturbations

We repeat the stability analysis by Cunningham (1976) and D’Alessio et al.
(1998), but using a form that will make generalization to a more complicated
form of $F_{\rm irr}$ straightforward. We perturb by setting $T\rightarrow
T+\delta T$, and similarly $F_{\rm irr}\rightarrow F_{\rm irr}+\delta F_{\rm
irr}$. Linearizing the thermal equation yields

$\displaystyle{{\partial}\over{\partial}t}{\delta T\over T}={1\over t_{\rm
th}}\left({\delta F_{\rm irr}\over{F_{\rm irr}}}-4{\delta T\over T}\right)\ .$
(14)

The power-law equilibrium profile, together with a constant $\chi$, yield

$\displaystyle{{\delta F_{\rm irr}}\over{F_{\rm irr}}}$ $\displaystyle=$
$\displaystyle\left({{d\over{d\ln r}}{H\over
r}}\right)^{-1}\,{{\partial\over{\partial\ln r}}\left({{\delta H}\over
H}{H\over r}\right)}$ (15) $\displaystyle=$ $\displaystyle{7\over
4}{\partial\over{\partial\ln r}}{{\delta T}\over T}+{1\over 2}{{\delta T}\over
T}\,,$

where the first term on the right-hand side bestows the travelling wave nature
for the perturbation, while the second term provides positive feedback. Now
the thermal equation becomes

$\displaystyle{{\partial}\over{\partial}t}{\delta T\over T}={1\over t_{\rm
th}}\left({7\over 4}{{\partial}\over{\partial}\ln r}{\delta T\over T}-{7\over
2}{\delta T\over T}\right)\,.$ (16)

This equation admits a decaying travelling wave solution. We consider a
complex $\delta T$ of the form

$\displaystyle{{\delta T}\over T}\propto e^{st+ik\ln r}\ ,$ (17)

with wavenumber $k$ (a real number, not restricted to integers) and growth-
rate $s$ (complex-valued). The physical perturbations are understood to be the
real parts of these. Inserting this into Equation (16) then yields the growth
rate

$\displaystyle s={1\over t_{\rm th}}\left({7\over 4}ik-{7\over 2}\right)\,,$
(18)

or, the thermal perturbation propagates inwards with a phase speed ${d\ln
r\over dt}|_{\rm phase}=-{1\over k}{\rm Im}[s]=-{7\over 4}{1\over t_{\rm
th}}$, and decays at a rate ${\rm Re}[s]=-{7\over 2}{1\over t_{\rm th}}$
(D’Alessio et al., 1998).

Why is the disk stable? Naively, one might imagine that a hotter region rises
up in scale height and can therefore intercept more stellar flux. However,
under the assumption that $\chi$ is constant, stellar heating scales with
temperature to the half power ($F_{\rm irr}\propto H\propto h\propto
T^{1/2}$), while cooling scales with a higher power ($\propto T^{4}$). So the
perturbation decays in time.

To understand why perturbations propagate inwards, we consider the disk
surface in the presence of a local positive temperature perturbation. The
inner half of the affected region has a steeper grazing angle than in
equilibrium, and so is heated more. This raises the local surface under
hydrostatic adjustment. Conversely, in the outer half, the surface drops. As a
result, the perturbation moves inwards.

##  5 Instability for Realistic Surface

A crucial effect not accounted for in the simplistic surface model is that the
depth of starlight penetration changes when the slope of the optical surface
varies (Fig. 1). That strengthens the heating response and leads to an
instability.

Here, we locate the disk surface using the definition that the optical depth
to the star is unity. The optical depth along a slanted ray from the star that
has an inclination angle $\theta_{H}$ is (Watanabe & Lin, 2008).

$\displaystyle\tau_{\rm slant}$ $\displaystyle=$
$\displaystyle\left(1+\tan^{2}\theta_{H}\right)^{1/2}\int_{0}^{r}\kappa_{V}\rho_{\rm
dust}(r^{\prime},z=\tan\theta_{H}r^{\prime})dr^{\prime}$ (19)
$\displaystyle\approx$ $\displaystyle\int_{0}^{r}\kappa_{V}\rho_{\rm
dust}(r^{\prime},z=\theta_{H}r^{\prime})dr^{\prime}\ ,$

where the approximate sign holds when $\theta_{H}\ll 1$. With the density
profile given by eq. (4), and the vertical optical depth
$\tau_{V}=\kappa_{V}\Sigma_{\rm dust}/2$, we have the following equation to
determine the more realistic optical surface,

$\displaystyle 1=\tau_{\rm
slant}\approx\int_{0}^{r}dr^{\prime}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}h(r^{\prime})}e^{-{1\over
2}{H(r)^{2}\over r^{2}}{r^{\prime 2}\over h(r^{\prime})^{2}}}\,.$ (20)

This is our revised equation for the optical surface. Given a disk opacity and
temperature profile, one can use it to solve for $H(r)$.

###  5.1 An approximate form and qualitative discussion

To better understand the implications of eq. (20), we derive a simpler form
that is valid for $\chi\equiv{H\over h}\gtrsim$ a few. The optical depth along
a given line of sight ($\theta_{H}$) is mostly produced by material close to
the optical surface, so we can write $\tau_{\rm slant}\approx\rho_{\rm
dust}(r,H)\,\kappa_{V}\,\Delta r$, where $\Delta r$ is the local density scale
height as experienced by the slanted ray. Expressing density as $\rho_{\rm
dust}(r,H)\propto e^{-H^{2}/(2h^{2})}\propto
e^{-r^{2}\theta_{H}^{2}/(2h^{2})}$ (eq. 4), we have

$\displaystyle{\Delta r}=\left.\left({d\ln\rho_{\rm dust}\over dr}\right|_{\rm
constant\,\theta_{H}}\right)^{-1}={r\over\chi^{2}\gamma}\ ,$ (21)

after defining a flaring index $\gamma$,

$\displaystyle\gamma\equiv{d\ln(h/r)\over d\ln r}\ .$ (22)

This index equals $2/7$ for the power-law disk (eq. 11). Setting $\tau_{\rm
slant}=1$ then yields our desired approximate form for the optical
surface,333Equation (23) may also be derived directly from equation (20) as
follows: assuming $\chi$ is large, the integrand in equation (20) is
exponentially suppressed unless $r^{\prime}$ is close to $r$. Therefore we may
approximate $r^{\prime}\approx r$ outside of the exponential, and inside of
the exponential we may set ${r^{\prime 2}\over
h(r^{\prime})^{2}}\approx{r^{2}\over h(r)^{2}}(1-2\gamma{r^{\prime}-r\over
r})$, in which case the integral yields equation (23) (Garaud & Lin, 2007).

$\displaystyle{\chi^{2}\over 2}e^{\chi^{2}\over
2}\approx{\tau_{V}\over\sqrt{2\pi}}{1\over\gamma\,h/r}\,.$ (23)

Although much simpler than eq. (20), this form still must be solved
numerically. However, it makes explicit the dependencies of $\chi$ on disk
properties such as the optical depth and the flaring angle, key physical
elements for the instability.

In the following, we discuss the main features of this model of a realistic
surface.

1\. Equilibrium disk. One can obtain the equilibrium disk by iteratively
solving the optical surface equation (eq. 20) and the steady-state thermal
equation. Appendix B.2 explains how we perform this procedure and achieve
convergence. The resulting value of $\chi$ for our fiducial disk is shown in
Fig. 2. It drops only modestly over a large stretch of the disk. This can be
easily understood from eq. (23): for $\chi\gg 1$, we can neglect terms outside
the exponential to find that $\chi$ depends logarithmically on $\tau_{V}$,
$\chi\sim\left(2\ln\tau_{\rm V}\right)^{1/2}$. Such a weak dependence on the
surface density arises due to the rapid fall-off of density with vertical
height, and allows CG97 to assume a constant $\chi$ throughout the disk and to
derive a power-law equilibrium solution (eq. 11).

2\. Instability. The most interesting implication of eq. (23), however, lies
in the relation between $\chi$ and the flaring index $\gamma$. When the
flaring index is larger, $\chi$ is smaller because the starlight shines closer
to the disk’s surface normal (larger grazing angle) and so can penetrate
deeper towards the midplane. It is this dependence on the grazing angle that
is crucial for the irradiation instability (§2).

To illustrate this point, we retain the key factors in eq. (23) to recast it
as

$\displaystyle e^{-{\chi^{2}\over 2}}\propto{{\partial}(h/r)\over{\partial}r}\
.$ (24)

A thermal perturbation of the spatial form $e^{ik\ln r}$ (with $k\gg 1$) then
perturbs the surface as444For a more exact form, c.f. eq. (B13).

$\displaystyle\delta\chi\approx-{{ik}\over{\chi\gamma}}{\delta h\over h}\ .$
(25)

This in turn affects the heating rate as (using eq. 2)

$\displaystyle\delta F_{\rm irr}$ $\displaystyle\propto$
$\displaystyle{{\partial}\over{\partial}\ln r}\delta
H={{\partial}\over{\partial}\ln r}\left(\chi\delta h+h\delta\chi\right)$ (26)
$\displaystyle=$ $\displaystyle ik\left(\chi\delta h+h\delta\chi\right)\,.$

While the first term in the brackets is present in the model of a simplistic
surface (constant $\chi$), the second term is not. It describes the change in
the penetration depth as the grazing angle varies. Inserting eq. (25) into
this expression, we find a positive definite contribution to the
heating,555For a more exact form, c.f. eq. (B14).

$\delta F_{\rm irr}\propto ik\chi\delta h+{{k^{2}}\over{\chi\gamma}}\delta
h\,.$ (27)

Compared with eq. (15), the second term is new. At sufficiently high
wavenumbers, it can overcome damping by radiative cooling and lead to a new
instability, the irradiation instability.

3\. The Smearing Length. There is a unique scale length associated with
irradiation. Thus far, we we have assumed that light is absorbed at a well
defined surface (the optical surface), But in truth, starlight is deposited
over a distance $\Delta r$ (eq. 21). We call the latter the ‘smearing length,’
and define an associated wavenumber

$\displaystyle k_{\rm smear}\equiv{r\over\Delta r}=\chi^{2}\gamma\ .$ (28)

With this new quantity, eq. (27) now reads

$\delta F_{\rm irr}\propto\left(i+{k\over{k_{\rm
smear}}}\right)\,k\chi\,\delta h\,.$ (29)

For $k\ll k_{\rm smear}$, the simplistic surface model prevails and thermal
perturbations lose out to radiative cooling. Modes with $k\sim k_{\rm smear}$
or larger, on the other hand, can be destabilized by the changes in the
penetration depth of the starlight.

Moreover, mode growth rates are affected by the smearing length. While the
above simple expression indicates that mode growth rate rises with $k$ as
$k^{2}$, a more careful analysis (Appendix B.4) that accounts for a finite
smearing length yields a saturated growth rate for modes with $k\gg k_{\rm
smear}$ (see text below).

###  5.2 Linear Perturbations

After the above qualitative discussions, we now present results from more
rigorous derivations. We call the following model ‘analytical,’ to distinguish
it from the RADMC simulations.

We employ equations (1), (2), (3), and (20) to study the stability of a
thermal perturbation. We first derive an approximate dispersion relation, then
compare it against exact numerical solutions. Our analysis shows that
perturbations of certain wavelengths indeed grow in amplitude, on a timescale
that is of order the local thermal time.

Assuming a space/time dependence of the form $e^{st+ik\ln r}$ (eq. 17), we
derive in Appendix B.4 an approximate analytic expression for the growth rate,

$\displaystyle s_{\rm grow}\equiv{\rm Re}[s]={1\over t_{\rm
th}}{k^{2}(\chi^{2}-8)-7k_{\rm smear}^{2}\over 2(k_{\rm smear}^{2}+k^{2})}\,,$
(30)

where $k_{\rm smear}=\chi^{2}\gamma$ (eq. 28). In the limit of long
wavelengths ($k\ll k_{\rm smear}$), variations of $\chi$ are relatively
insignificant, and so one recovers that the wave damps at the rate $s_{\rm
grow}=-7/(2t_{\rm th})$ (eq. 18, and D’Alessio et al., 1998).

Equation (30) allows unstable modes whenever $\chi>\sqrt{8}$, i.e., when the
disk is sufficiently opaque that the optical surface lies well above the gas
scale height. As Fig. 2 shows, this condition is satisfied throughout our
fiducial disk. The unstable waves have short wavelengths,

$k\geq k_{\rm min}\equiv\sqrt{{7\over{\chi^{2}-8}}}k_{\rm smear}\,,$ (31)

with the value of $k_{\rm min}$ ranging from $4.3$ to $4.8$ in our fiducial
disk (see Fig. 5). For reference, the wavelength of a $k=2.7$ wave spans one
decade in radius. So these unstable waves have $\sim 2$ or more wavelengths
per radial decade.

The growth rates of unstable waves first rise with $k$ as $k^{2}$, before
saturating to a constant value for $k\gg k_{\rm smear}$. This saturation is
related to the smearing length, i.e., the finite spread of the optical
surface. When the wavelength is much shorter than this length, changes in
stellar heating are subdued relative to that for a razor-sharp optical
surface. As a result, the growth rate saturates.

Figure 3:  Evolution of linear perturbations, of the form $e^{ik\ln r}$, in
the model with a realistic optical surface. The top panel is for a $k=2$
perturbation, while the bottom is for $k=6$. Snapshots are taken every $1/4$
of the thermal time at 10AU, with the initial ones being the lightest in
color. The $k=2$ perturbation decays in time, while the $k=6$ one grows.
Figure 4: Growth rates for linear perturbations of the form $e^{ik\ln r}$, in
the model with a realistic surface. The growth rates are in units of inverse
local thermal time, with different colors standing for different $k$-values.
The thick curves show the exact numerical result, while the thin dotted curves
present results of our analytic approximation (eq. 30). Throughout our
fiducial disk, $\chi\geq\sqrt{8}$, and there are unstable modes.

To confirm these analytical findings, we obtain mode growth rates by
numerically integrating the relevant equations. The numerical details are in
Appendix B.1. Fig. 3 shows the results of two such integrations. Starting from
an equilibrium disk (Appendix B.2), we impose small initial sinusoidal
perturbations of the form $e^{ik\ln r}$. We observe that the first
perturbation (with $k=2$) damps with time as the wave travels inward, while
the second one ($k=6$) grows.

We can extract growth rates in integrations like these, following the
procedure in Appendix B.3. The results for some low $k$ waves are plotted in
Fig. 4. We find that the analytical expression (eq. 30) agrees with the
numerical results qualitatively. Both suggest that $k\geq 4$ waves should be
unstable over much of our fiducial disk, with a growth time that is of order
the local thermal time.

##  6 Instability also found by RADMC

We turn to the radiative transfer code RADMC. We do so for multiple purposes.
One is to substantiate our analytical results. Another is to use RADMC to
remedy a major shortcoming in the analytical approach, i.e., horizontal
averaging (see below). An additional advantage of RADMC is that it allows us
to go beyond the two-frequency approximation.

Given an assumed dust density field, RADMC uses the Monte Carlo method to
follow photon absorption and emission, and to determine the local equilibrium
temperature by balancing energy gain and loss. In principle, it can be
combined with the equation of hydrostatic equilibrium to obtain, iteratively,
the equilibrium disk profile. This has been attempted by a number of studies
(e.g., Dullemond & Dominik, 2004b; Min et al., 2009; Siebenmorgen & Heymann,
2012; Ueda et al., 2019), either with RADMC or an analogous code. But the
procedure does not converge for optically thick disks. We show in Appendix C
that this non-convergence is partly caused by a numerical (non-physical)
instability. So non-convergence does not prove that the disk is truly
unstable. To tease out the the physical instability requires some finesse.

###  6.1 Simulation Setup

We adopt the same fiducial disk as before. The density profile is assumed to
stay in hydrostatic equilibrium (eq. 4) as the midplane temperature $T=T(r,t)$
evolves. RADMC treats the radiation field at multiple wavelengths. The opacity
law is set following eq. (7).

To trace the thermal evolution, we replace the stellar heating term in eq. (1)
by

$\displaystyle{1\over 2}F_{\rm irr}[T]\rightarrow\sigma_{\rm SB}\left(T_{\rm
RADMC}[T]\right)^{4}\,,$ (32)

where $T_{\rm RADMC}[T]$ is the midplane temperature profile output by
RADMC.666We use the midplane temperature, because in RADMC simulations, the
disk is vertically isothermal below the optical surface, and the midplane
temperature reflects the amount of stellar heating for the disk. In contrast,
the optical surface is hotter (Chiang & Goldreich, 1997). We use this as a
proxy for stellar heating. The thermal equation becomes

$\displaystyle{3\over 8}\tau_{\rm mm}c_{p}\Sigma_{\rm
gas}{{\partial}T\over{\partial}t}=\sigma_{\rm SB}\left(\left(T_{\rm
RADMC}[T]\right)^{4}-T^{4}\right)\,.$ (33)

At equilibrium, $T=T_{\rm RADMC}[T]$, as is desired. Whereas previous works
have ignored the left-hand term when iteratively searching for equilibrium
solution, we integrate the full equation forward in time. At each time-step,
we inject 20 million stellar photon packets from the origin. This provides a
sufficiently accurate map of the illumination pattern. We choose a time-step
that is a small fraction of the thermal time (typically 4% of the thermal time
at 10AU).

Before we proceed to present results of these integrations, we comment on
three major issues. The first two cast some doubts on the RADMC results, the
third relates to a major improvement of RADMC over our analytical study.

First, the inner boundary of our disk is set to be 1AU. Between 1 and 1.3AU,
we freeze the disk profile to avoid a puffed-up inner rim, which would
otherwise cast a shadow further out. Although such a shadow may indeed be
realistic, we do not wish its effect to pollute what happens in the simulation
domain. Second, eq. (33) itself contains an important shortcoming: excessive
thermal diffusion. In realistic situations, temperature gradients in disks
should be communicated on diffusive timescales. However, the way we hijack
RADMC for our purpose effectively assumes that both vertical and horizontal
heat transfer are instantaneous. The former shortcoming (neglecting vertical
diffusion) produces a greater perturbation on the midplane temperature than
reality, boosting the instability; the latter one (horizontal diffusion), on
the other hand, leads to an excessive damping of the instability. Clearly, a
more thorough study that considers time-dependent radiative transfer is
required.

The third major issue is what we term ’horizontal averaging’. This is
different from the horizontal diffusion discussed in the last paragraph. As
grains at the perturbed optical surface receive more (less) stellar
irradiation, they light up as brighter (dimmer) bulbs for the disk down below.
In our analytical study, we have assumed that the midplane disk can only see
bulbs lying directly above it (heating is only affected by the local gradient,
eq. 2). This assumption fails if the distance between the bulbs and the
midplane (more accurately, the vertical photosphere in thermal radiation) is
greater than the radial wavelength, $H\geq r/k$. In this case, the midplane
feels a reduced heating (cooling) due to horizontal averaging of the different
light bulbs, and we expect the instability to be quenched. This effect was
modelled in Watanabe & Lin (2008) by manually averaging the perturbations over
a radial distance; in RADMC simulations, it is automatically captured.

To order of magnitude, horizontal averaging limits unstable waves to those
with

$k\leq k_{\rm max}\equiv{r\over H}\,.$ (34)

Together with the lower limit on unstable $k$ (eq. 31), this defines the
region in disks where irradiation instability can occur. For our fiducial
disk, we expect unstable waves only inward of $\sim 30$AU (Fig. 5).

Figure 5: The expected unstable region for our fiducial disk. Wavenumbers of
unstable modes have to fall within the shaded zone: $k_{\rm min}\leq k\leq
k_{\rm max}$ (eq. 31 & eq. 34). This occurs inward of $\sim 30$AU. For
reference, $k=2.7$ corresponds to 1 wavelength per radius decade.

###  6.2 Linear Evolution

Figure 6:  Similar to Fig. 3 but now obtained using RADMC integrations. We
observe that the $k=2$ perturbation damps, while the $k=6$ perturbation grows,
as is the case in our analytical model. However, growth for the latter is now
restricted to only the inner disk ($r<30$AU).

We impose a small sinusoidal perturbation of the form $e^{ik\ln r}$, on the
equilibrium disk (cf. §5). The initial linear evolution is shown in Fig. 6.
Comparing with the analogous plot for our analytical model (Fig. 3), we see
that the $k=2$ perturbation decays in both calculations. On the other hand,
while our analytical model predicts growth for the $k=6$ perturbation at all
radii, RADMC runs show that it only experiences growth inside 30AU. This
reflects the short-coming of the analytical model: it does not account for
horizontal averaging, a problem that is more severe in the outer disk.

Figure 7:  Growth rates of linear perturbations extracted from RADMC
simulations (jagged thick curves), and from the analytical model (smooth thick
curves, from Fig. 4). The former show in general lower growth rates – all
modes are stable outside 30AU – due to both horizontal averaging and
horizontal diffusion (see text). RADMC results inward of $\sim 1.3$AU, where
we freeze the disk profile, should be ignored.

We design a special numerical procedure to extract mode growth rates from
RADMC simulations. This is detailed in Appendix B.3. The results are shown in
Fig. 7 for a few wavenumbers. In this exercise, in order to use RADMC to
accurately track changes in heating associated with a small perturbation, we
launch a very large number of photon packets (20 billion). Despite this, the
results still appear somewhat jagged.

We are able to confirm that RADMC yields qualitatively similar growth rates as
those from the analytical model (thick lines in Fig. 4, also reproduced here).
But there are two notable differences. First, RADMC shows that the outer disk
(beyond 30AU) is stable to perturbations. As we explain above, this is related
to the issue of horizontal averaging in the outer region. Second, the growth
rates for unstable modes are somewhat lower in RADMC, indicating enhanced
damping. We believe this may be related to the fact that the RADMC procedure
assumes instantaneous horizontal diffusion. A fully time-dependent radiative
transfer code, together with a hydrodynamic solver, will be necessary to
provide more accurate answers.

###  6.3 Nonlinear Development

Figure 8: Late-time evolution of the midplane temperature (top panel) and the
optical surface (bottom panel), obtained from RADMC simulations. Starting from
the equilibrium profile (labelled $t=0$), the increasingly dark curves display
snapshots at 5, 6, and 7 thermal times at $10$ AU. The disk inward of 1.3AU is
frozen in the simulations. One observes instability only inward of $\sim 30$
AU.

RADMC also allows us to follow the instability to its nonlinear stage. For the
following runs, we use 20M photon packets, and have checked to confirm that
the behavior does not depend on numerical parameters such as timestep, grid
spacing, initial conditions, and inner/outer boundary conditions.

Starting from an equilibrium profile, small perturbations (seeded only by
numerical noise) grow in time. After around 5 thermal times at 10AU, the
system reaches a quasi-steady-state, in which waves are generated around
$30$AU, and grow in amplitude as they propagate inwards. A few snapshots are
presented in Fig. 8. Globally, the most prominent waves have wavenumbers
$k\sim$ 4-6 (corresponding to about two wavelengths per radius decade), as
expected from their growth rates and horizontal averaging (Fig. 5) . While the
disk outside $30$AU remains largely unperturbed, inside $30$AU these waves
cause order-unity perturbations in the midplane temperature. The bottom panel
of Fig. 8 shows the optical surface at $0.5\mu$m. This looks like a staircase:
the front edge of each stair is illuminated by the star (“stair riser”), and
behind that the disk falls into shadow and the optical surface appears flat in
$H/r$ (“stair tread”).

In the nonlinear stage, the temperature waves are asymmetric with steeper sun-
facing edges, as is also found by Watanabe & Lin (2008). This may arise
because at large amplitudes, the local flaring index ($\gamma$) of the front
edge is higher, allowing waves of higher wavenumber to be unstable.777The
inner edge of a wave could also be thought of as a puffed-up inner rim,
casting a long shadow on the disk behind. As such, perturbations with
initially long wavelengths can acquire sharp spatial gradients, leading to
large pressure variations across radial scales as small as the disk scale
height. In fact, as our RADMC procedure introduces enhanced horizontal heat
diffusion, real disks may harbour waves that are even larger in amplitudes and
sharper in scale.

##  7 Re-cap of Technical Results

We summarize results obtained in previous sections.

We found that passive disks are susceptible to an irradiation instability. A
slight thermal perturbation in such a disk can lead to the disk receiving even
more stellar heating, thereby initiating unstable growth. We examined this
instability using three models of increasing complexity and realism.

In the simplest model, we assumed that the height of the optical surface (the
disk layer at which stellar photons are absorbed) varies in proportion to the
gas scale height when the disk is thermally perturbed. We found that thermal
perturbations propagate inwards and decay as they do so, in agreement with
D’Alessio et al. (1999). The modes are damped because the increase in stellar
insolation in this model is insufficient to counteract the damping due to
blackbody cooling.

In contrast, in a model where the optical surface is self-consistently
determined, we found that sinusoidal perturbations can grow as they propagate
inward. The key physics is that as a thermal perturbation increases the scale
height in a disk annulus, the inner half of the annulus acquires a steeper
flaring, and the outer half a more gentle one. Starlight now shines more
overhead for the inner half and so can penetrate into deeper layers;
conversely, it is absorbed more shallowly in the outer part. These geometric
changes allow the annulus to intercept more starlight, more so than in the
first model, and give rise to instability. Thermal perturbations grow and
travel inwards in of order the local thermal time.

In terms of unstable wavelength, very long waves are stable because they cause
little change in the surface slope; only waves with wavelengths comparable or
shorter than the so-called ‘smearing length’, the slant length over which the
star deposits its energy, can grow. More optically thick disks have shorter
smearing lengths and therefore harbour a broader spectrum of unstable waves.
For our fiducial disk, the unstable waves have wavenumber $k$ $\geq 4$, i.e.,
with $\sim 2$ or more wavelengths per decade in radius.

To confirm those analytical results, we retooled the radiative transfer code
RADMC to track the stellar heating. This also allowed us to address the issue
of horizontal averaging, the complication that arises when the disk heating is
not solely determined by the slope of the local optical surface, but also by
nearby hot and cold patches on the optical surface. Using both scaling
arguments and RADMC simulations, we found unstable modes exist only inside
$\sim 30$AU, or, where the disk is geometrically thin ($H\leq r/k$).

We also used RADMC to study the nonlinear evolution of the irradiation
instability. Within a few thermal times, waves generated far from the star
have propagated to the inner region and have grown to order unity amplitudes.
The front edges of these waves can be as sharp as a pressure scale height,
suggesting that they may be effective in trapping dust grains.

###  7.1 Relation to Watanabe & Lin

The irradiation instability was first discussed in Watanabe & Lin (2008).
Calculating the optical surface self-consistently (as in our ‘realistic
surface’ model), they numerically integrated the thermal equation and observed
that thermal waves are driven to large amplitudes as they travel inwards. But
unlike our RADMC simulations here, they found instability throughout the whole
disk. This difference may partly be explained by the different ways we account
for horizontal averaging: while RADMC naturally captures this phenomenon,
Watanabe & Lin (2008) performed a radial average of the surface re-radiation.
Moreover, compared to their numerical results, our analysis here has the
advantage of elucidating the origin of the instability, and identifying the
unstable wavenumbers.

##  8 Examining Assumptions

There remain a wide range of uncertainties associated with the irradiation
instability. In this preliminary exploration, we have had to make a string of
assumptions. Here, we discuss some of these.888Also see Appendix D which
validates the assumption of rapid gas-dust thermal coupling. We suggest that
modifying them can lead to an array of interesting results. But the essence of
the irradiation instability, we feel, should prevail.

###  8.1 Hydrodynamics

We simplify the hydrodynamical response of the disk by assuming that the scale
height reacts instantaneously to the midplane temperature, the perturbation is
axisymmetric, and the surface density does not evolve. These are questionable
for multiple reasons.

First, they are predicated on the assumption that the orbital time is shorter
than the thermal time, which is only true inside of $\sim 50$AU in our
fiducial disk (Fig. 2). What happens at larger radii is uncertain.
Perturbation analyses by Dullemond (2000) and Chiang (2000) show that a
different instability may operates in that limit. But a more careful analysis,
including effects such as starlight penetration depth, smearing length, etc.
is needed. In addition, perturbations are likely no longer axisymmetric when
the orbital time is long.

Second, even in the limit of long thermal time, hydrodynamical effects can
dramatically alter the story. A local heating will not only change the gas
scale height, but can also expel material to neighbouring annuli. This may set
up meridional flows that advect heat and dust. Vortices or turbulence may also
ensue (via, e.g., the Rossby wave instability, Lovelace et al., 1999;
Lithwick, 2009). Investigating these effects, both in terms of how they impact
the irradiation instability and of how they impact the long-term disk
evolution, requires numerical simulations that co-evolve hydrodynamics and
radiation. We are currently pursuing this course. Another motivation for a
full hydro+radiation treatment is to improve on our treatment of thermal
diffusion. Thermal perturbations in optically thick disks should propagate as
diffusive waves. However, even in our most sophisticated account of the
thermal physics (the RADMC approach), this is not accurately captured. Our
retooling of RADMC accelerates horizontal diffusion artificially. We suspect
that this leads to enhanced damping for the waves, as well as reduced
amplitudes during the nonlinear evolution. We have also ignored vertical
thermal waves, which may have the potential to weaken the instability (but see
Appendix A.3).

###  8.2 Dust Movement

Figure 9: An initial foray into dust settling. The set-up is identical to that
in Fig. 8 except that we imitate dust settling by removing all dust grains
above $2$ scale heights. The instability now extends to larger distances (out
to $\sim 100$AU).

We assume that dust and gas are co-spatial, with a constant mass opacity
everywhere in the disk. This is perhaps our most problematic assumption. It is
expected that dust grains evolve in at least three dimensions: they settle
vertically, migrate radially, and their sizes can change due to conglomeration
and fragmentation. This can change the opacity relative to what we adopt here.

Dust opacity, especially that in the optical, is the deciding factor for the
irradiation instability. The discussion following eq. (30) indicates that
instability requires a fairly optically thick disk. Opacity in the optical is
mostly contributed by micron-sized grains. These grains are known to settle
quickly from these heights (Dullemond & Dominik, 2004b), though even weak
turbulence or circulation can disrupt this process. While the jury is still
out on the vertical distribution of these small grains, we explore a scenario
where settling has occurred. In Fig. 9, we repeat our RADMC calculations but
now with all dust above $2h$ removed. The optical surface now lies much closer
to the midplane ($H\leq 2h$). This mitigates the negative impact of horizontal
averaging (eq. 34). And our RADMC simulations show that the instability now
extends to almost the entire disk (out to $\sim 100$AU).

In summary, it appears that a complete picture of the irradiation instability
can only emerge after we understand dust physics, an endeavour that is further
complicated by the fact that the hydrodynamical response of the disk affects
dust settling and fragmentation (§8.1).

###  8.3 Inner Region

We simplified the physics in the inner region substantially, by ignoring the
presence of a hot rim and by assuming the star is point-like. Moreover, we
have ignored viscous heating, which may dominate over starlight heating in the
inner region. Here we briefly discuss how these three effects may impact the
irradiation instability.

A hot rim is expected to form where the inner disk sees the star unobstructed
(Natta et al., 2001; Dullemond et al., 2001). Observed around T Tauri stars
(Muzerolle et al., 2003; Akeson et al., 2005) and Herbig Ae/Be stars (Natta et
al., 2001), such a rim casts a long shadow on the disk behind (Dullemond et
al., 2001; Vinković et al., 2006; Siebenmorgen & Heymann, 2012; Flock et al.,
2016). In this report, we have opted to simplify the picture by freezing the
disk profile near the inner edge. But our (additional) explorations with RADMC
often find that disks near the shadow terminator are unstable, in agreement
with some previous work (Siebenmorgen & Heymann, 2012; Flock et al., 2016;
Ueda et al., 2019). We suspect that it is related to the irradiation
instability, but have yet to provide firm evidence.

We also simplified the star into a point-like light source. This is adequate
as long as $0.4R_{*}\ll{2\over 7}{H}$ (Kusaka et al., 1970; Chiang &
Goldreich, 1997), where $R_{*}$ is the stellar radius, or outside $0.4$AU for
a star with a radius of $2.5R_{\odot}$. Inside this region, we expect the
irradiation instability to be suppressed by the finite-source effect, and this
is indeed demonstrated by the numerical experiments of Watanabe & Lin (2008).

Viscous heating associated with mass accretion can compete against stellar
heating only in the very inner disk (inward of an AU). The irradiation
instability should be suppressed there (Watanabe & Lin, 2008).

##  9 Conclusions

We have developed a model to study the irradiation instability in passively
heated disks. Our approach, combining analytical and numerical tools, reveals
the origin of the instability, identifies the unstable wavelengths and obtains
relevant growth-rates. Our fiducial disk is unstable inside of $30$AU, and
produces large-amplitude inwardly-propagating thermal waves. A preliminary
exploration shows that the instability can extend to much further reaches, if
some degree of dust settling is included.

Future work should relax many of the approximations we have made here. Of
particular urgency are the assumptions of instantaneous hydrostatic
adjustment, vertical isothermal structure, and dust evolution.

If the instability proves robust, it has many important potential
implications. For example, the large amplitude waves could be the forebears
for the gaps and rings observed by ALMA. These latter features are observed at
radial distances of 10-100AU, and have radial scales that are compatible with
the unstable wavelengths we obtain here ($\sim 2$ wavelengths per decade in
radius). The narrowness of the bright rings (Dullemond et al., 2018) is also
consistent with the nonlinear development of the irradiation instability.
Moreover, these waves can affect dust migration and conglomeration, which in
turn feeds back on the instability. Last but not least, the same waves could
also be hydrodynamically unstable, and produce vertical circulation, vortices
and turbulence. This could provide the much needed viscosity for disk
accretion. Clearly, much exciting work lies ahead.

Y.W. acknowledges support from NSERC. Y.L. acknowledges NASA grant NNX14AD21G
and NSF grant AST-1352369. We thank Kees Dullemond for making the RADMC-3D
code public and highly usable, and an anonymous referee for a careful and
insightful review.

## References

  * Akeson et al. (2005) Akeson, R. L., Walker, C. H., Wood, K., et al. 2005, ApJ, 622, 440, doi: 10.1086/427770
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andrews (2020) Andrews, S. M. 2020, arXiv e-prints, arXiv:2001.05007.  https://arxiv.org/abs/2001.05007
  * Andrews & Birnstiel (2018) Andrews, S. M., & Birnstiel, T. 2018, Dust Evolution in Protoplanetary Disks, ed. H. J. Deeg & J. A. Belmonte (Springer), 136, doi: 10.1007/978-3-319-55333-7_136
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
  * Bae et al. (2017) Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201, doi: 10.3847/1538-4357/aa9705
  * Bai (2014) Bai, X.-N. 2014, ApJ, 791, 137, doi: 10.1088/0004-637X/791/2/137
  * Bai (2015) —. 2015, ApJ, 798, 84, doi: 10.1088/0004-637X/798/2/84
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76, doi: 10.1088/0004-637X/769/1/76
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Baruteau et al. (2014) Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 667, doi: 10.2458/azu_uapress_9780816531240-ch029
  * Bowler (2016) Bowler, B. P. 2016, Publications of the Astronomical Society of the Pacific, 128, 102001, doi: 10.1088/1538-3873/128/968/102001
  * Burrows et al. (1996) Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437, doi: 10.1086/178156
  * Chiang (2000) Chiang, E. 2000, PhD thesis, CALIFORNIA INSTITUTE OF TECHNOLOGY
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Chiang et al. (2001) Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., et al. 2001, ApJ, 547, 1077, doi: 10.1086/318427
  * Cleeves et al. (2016) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110, doi: 10.3847/0004-637X/832/2/110
  * Cunningham (1976) Cunningham, C. 1976, ApJ, 208, 534, doi: 10.1086/154636
  * D’Alessio et al. (2006) D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314, doi: 10.1086/498861
  * D’Alessio et al. (1998) D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411, doi: 10.1086/305702
  * D’Alessio et al. (1999) D’Alessio, P., Cantó, J., Hartmann, L., Calvet, N., & Lizano, S. 1999, ApJ, 511, 896, doi: 10.1086/306704
  * Dipierro et al. (2016) Dipierro, G., Laibe, G., Price, D. J., & Lodato, G. 2016, MNRAS, 459, L1, doi: 10.1093/mnrasl/slw032
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127, doi: 10.3847/1538-4357/aa72f2
  * Dong et al. (2015) Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93, doi: 10.1088/0004-637X/809/1/93
  * Dullemond (2000) Dullemond, C. P. 2000, A&A, 361, L17
  * Dullemond & Dominik (2004a) Dullemond, C. P., & Dominik, C. 2004a, A&A, 421, 1075, doi: 10.1051/0004-6361:20040284
  * Dullemond & Dominik (2004b) —. 2004b, A&A, 417, 159, doi: 10.1051/0004-6361:20031768
  * Dullemond & Dominik (2005) —. 2005, A&A, 434, 971, doi: 10.1051/0004-6361:20042080
  * Dullemond et al. (2001) Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957, doi: 10.1086/323057
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool.  http://ascl.net/1202.015
  * Dullemond & Penzlin (2018) Dullemond, C. P., & Penzlin, A. B. T. 2018, A&A, 609, A50, doi: 10.1051/0004-6361/201731878
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46, doi: 10.3847/2041-8213/aaf742
  * Flaherty et al. (2018) Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117, doi: 10.3847/1538-4357/aab615
  * Flock et al. (2016) Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144, doi: 10.3847/0004-637X/827/2/144
  * Flock et al. (2015) Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68, doi: 10.1051/0004-6361/201424693
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155, doi: 10.3847/1538-4357/ab9641
  * Furlan et al. (2006) Furlan, E., Hartmann, L., Calvet, N., et al. 2006, ApJS, 165, 568, doi: 10.1086/505468
  * Garaud & Lin (2007) Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606, doi: 10.1086/509041
  * Gaudi (2021) Gaudi, B. S. 2021, arXiv e-prints, arXiv:2102.01715.  https://arxiv.org/abs/2102.01715
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Kenyon & Hartmann (1987) Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714, doi: 10.1086/165866
  * Klahr et al. (2018) Klahr, H., Pfeil, T., & Schreiber, A. 2018, Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation, ed. H. J. Deeg & J. A. Belmonte (Springer International Publishing AG), 138, doi: 10.1007/978-3-319-55333-7_138
  * Kusaka et al. (1970) Kusaka, T., Nakano, T., & Hayashi, C. 1970, Progress of Theoretical Physics, 44, 1580, doi: 10.1143/PTP.44.1580
  * Li et al. (2021) Li, J., Dempsey, A. M., Li, H., & Li, S. 2021, arXiv e-prints, arXiv:2102.02216.  https://arxiv.org/abs/2102.02216
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17, doi: 10.1088/0004-637X/811/1/17
  * Lithwick (2009) Lithwick, Y. 2009, ApJ, 693, 85, doi: 10.1088/0004-637X/693/1/85
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, Monthly Notices of the Royal Astronomical Society, 168, 603, doi: 10.1093/mnras/168.3.603
  * Min et al. (2009) Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155, doi: 10.1051/0004-6361/200811470
  * Muzerolle et al. (2003) Muzerolle, J., Calvet, N., Hartmann, L., & D’Alessio, P. 2003, ApJ, 597, L149, doi: 10.1086/379921
  * Natta et al. (2001) Natta, A., Prusti, T., Neri, R., et al. 2001, A&A, 371, 186, doi: 10.1051/0004-6361:20010334
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13, doi: 10.3847/1538-3881/ab16e9
  * Okuzumi et al. (2016) Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82, doi: 10.3847/0004-637X/821/2/82
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, The Astrophysical Journal, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Riols et al. (2020) Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95, doi: 10.1051/0004-6361/201937418
  * Safronov (1962) Safronov, V. S. 1962, Soviet Ast., 6, 217
  * Siebenmorgen & Heymann (2012) Siebenmorgen, R., & Heymann, F. 2012, A&A, 539, A20, doi: 10.1051/0004-6361/201118493
  * Simon et al. (2015) Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117, doi: 10.1093/mnras/stv2070
  * Suzuki et al. (2016) Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145, doi: 10.3847/1538-4357/833/2/145
  * Takahashi & Inutsuka (2014) Takahashi, S. Z., & Inutsuka, S. 2014, The Astrophysical Journal, 794, 55, doi: 10.1088/0004-637x/794/1/55
  * Ueda et al. (2021) Ueda, T., Flock, M., & Birnstiel, T. 2021, ApJ, 914, L38, doi: 10.3847/2041-8213/ac0631
  * Ueda et al. (2019) Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10, doi: 10.3847/1538-4357/aaf3a1
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132, doi: 10.3847/1538-4357/aa5d68
  * Vinković et al. (2006) Vinković, D., Ivezić, Ž., Jurkić, T., & Elitzur, M. 2006, ApJ, 636, 348, doi: 10.1086/497895
  * Wang & Goodman (2017) Wang, L., & Goodman, J. 2017, ApJ, 847, 11, doi: 10.3847/1538-4357/aa8726
  * Watanabe & Lin (2008) Watanabe, S., & Lin, D. N. C. 2008, ApJ, 672, 1183, doi: 10.1086/523347
  * Watanabe et al. (1990) Watanabe, S.-I., Nakagawa, Y., & Nakazawa, K. 1990, ApJ, 358, 282, doi: 10.1086/168985
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57, doi: 10.1093/mnras/180.2.57
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103, doi: 10.1051/0004-6361/201526538
  * Woitke et al. (2019) Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301, doi: 10.1088/1538-3873/aaf4e5
  * Wünsch et al. (2005) Wünsch, R., Klahr, H., & Różyczka, M. 2005, MNRAS, 362, 361, doi: 10.1111/j.1365-2966.2005.09319.x
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744

##  Appendix A Derivation of the Thermal Equation

###  A.1 Flux Splitting

In a passively irradiated disk, there are three temperatures of concern. The
first is that at the optical surface, $T_{s}$. This layer receives the stellar
flux $F_{\rm irr}$, then emits half down, half up. It has negligible thermal
inertia and simply acts as a beam splitter. We have

$\displaystyle F_{\rm irr}=2\sigma_{\rm SB}T_{s}^{4}\,.$ (A1)

The downward flux, $\sigma_{\rm SB}T_{s}^{4}$, is rapidly reprocessed at the
disk thermal photosphere which has a temperature $T_{\rm ph}$. This layer
emits a blackbody flux upwards, and transmits a flux $F_{\rm in}$ downwards
toward the optically thick midplane, leading to

$\sigma_{\rm SB}T_{s}^{4}=\sigma_{\rm SB}T_{\rm ph}^{4}+F_{\rm in}\,.$ (A2)

In steady state, $F_{\rm in}=0$ and the disk is vertically isothermal, with
midplane temperature $T_{c}=T_{\rm ph}=T_{s}=(F_{\rm irr}/2\sigma_{\rm
SB})^{1/4}$.

###  A.2 Simple Thermal Equation

We adopt a simple model for describing the thermal evolution of the disk. We
assume the downward diffusive flux determined by the equation of radiative
diffusion,

$\displaystyle F_{\rm in}={4\sigma_{\rm SB}\over 3}{{T_{\rm
ph}^{4}-T_{c}^{4}}\over\tau_{\rm mm}}\,,$ (A3)

where $T_{c}$ is the midplane temperature and $\tau_{\rm mm}$ is the optical
depth to the midplane for thermal radiation. We then assume that this inward
flux heats up gas near the midplane at the rate

$c_{p}{\Sigma_{\rm gas}\over 2}\,{{\partial T_{c}}\over{\partial t}}=F_{\rm
in}\,,$ (A4)

where $T_{c}$ is the midplane temperature and $\Sigma_{\rm gas}$/2 is the
surface density from the midplane outwards.

Here, we have adopted $c_{p}$ as the gas heat capacity. Watanabe et al. (1990)
have argued for an extra factor of $(\Gamma+1)/2\Gamma$, where $\Gamma$ is the
adiabatic index, when one also accounts for the vertical expansion of the disk
under stellar gravity. Moreover, when the disk’s thermal state is changed, its
rotational speed also has to adapt in order to stay in hydrostatic
equilibrium. In this exploration, we have ignored these corrections as they
only change the timescales by an order-unity factor.

Combining Equations (A1)–(A4), we land at our desired thermal equation

${1\over 2}(1+{3\over 4}\tau_{\rm mm})c_{p}\Sigma_{\rm gas}\,{{\partial
T_{c}}\over{\partial t}}={1\over 2}F_{\rm irr}-\sigma_{\rm SB}T_{c}^{4}\,.$
(A5)

In the main text, we drop the 1 in the brackets, as is valid for $\tau_{\rm
mm}\gtrsim 1$. This derivation ignores heat diffusion in the radial direction,
which is reasonable since the radial wavelength is typically much longer than
the scale height. It also does not properly account for the vertical
temperature profile, which we discuss next.

###  A.3 Vertical Thermal Waves

Equations (A3)– (A4) approximate the vertical transport of heat when the
timescale of variation is longer than the thermal time. Since unstable modes
have those two timescales comparable to each other, we examine here what
happens in the opposite limit. Our principal conclusions are (i) that vertical
thermal waves appear; and (ii) although those waves could potentially kill off
instability, they probably do not, particularly if the mm-sized dust has
settled.

We solve for the vertical temperature structure when the perturbation $\delta
T_{\rm ph}\exp(-i\omega t)$. Assume the disk is optically thick and energy
transport is by radiative diffusion. Energy conservation gives at any height

$\rho_{\rm gas}c_{p}{{\partial T}\over{\partial
t}}=\nabla\cdot\left({{16\sigma_{\rm SB}T^{3}}\over{3\kappa\rho_{\rm
dust}}}\nabla T\right)\,.$ (A6)

Here, gas provides the thermal capacity ($\rho_{\rm gas}$$c_{p}$ on the left-
hand side), while dust provides the opacity ($\kappa\rho_{\rm dust}$ on the
right-hand side).

It is preferable to use optical depth in lieu of the spatial coordinate,
$\tau_{z}=\int_{z}^{\infty}\kappa\rho_{\rm dust}dz$. Eq. (A6) now becomes

${{{3c_{p}}\over{16\sigma_{\rm SB}\kappa}}}{\rho_{\rm gas}\over\rho_{\rm
dust}}{{\partial T}\over{\partial
t}}={\partial\over{\partial\tau_{z}}}\left(T^{3}{{\partial
T}\over{\partial\tau_{z}}}\right)\,.$ (A7)

This is valid even if the dust has settled relative to the gas. But presently
we assume no settling. Note that eq. (A4) is essentially the the vertical
integral of eq. (A7), together with the approximation that
$T^{4}\propto\tau_{z}$ in the unperturbed state.

The solution for the temperature perturbation, subject to the boundary
conditions at the photosphere and at depth (i.e., $\delta T\rightarrow 0$ at
large $\tau_{z}$) is

$\delta T(\tau_{z},t)=\delta T_{\rm
ph}\exp\left(-{{\tau_{z}}\over{\tau^{\prime}}}\right)\cos\left(\omega
t-{{\tau_{z}}\over{\tau^{\prime}}}\right)\,,$ (A8)

where

$\displaystyle\tau^{\prime}=\left({32\sigma_{\rm SB}\kappa T^{3}\over
3c_{p}\omega}{\rho_{\rm dust}\over\rho_{\rm
gas}}\right)^{1/2}=\sqrt{8}{\tau_{z}(z=0)\over\sqrt{\omega t_{\rm th}}}\,.$
(A9)

Equation (A8) describes a wavetrain that propagates downwards. The amplitude
of the wavetrain decays exponentially. The key dimensionless parameter is

$\displaystyle{\tau_{z}\over\tau^{\prime}}\big{|}_{\rm
midplane}={1\over\sqrt{8}}\sqrt{\omega t_{\rm th}}\approx\sqrt{{7k\over
32}}=\sqrt{k/4.6}\ ,$ (A10)

where in the approximate equality we used the dispersion relation for
irradiated thermal waves (eq. 18). If that parameter is $\gtrsim{\pi\over 2}$,
then the vertical temperature profile will be oscillatory, which will likely
suppress the irradiation instability. Our unstable modes have $k\sim 6$, which
suggests that the instability is not suppressed. But clearly a more careful
treatment is needed to obtain more accurate values for the order-unity
coefficents.

One should also account for dust settling, because the mm-sized grains
responsible for the opacity tend to settle towards the midplane. The waves in
eq. (A8) are in $\tau_{z}$ rather than $z$. So if the grains settle into a
thin layer, then even if the parameter in eq. (A10) exceeds unity, the
temperature throughout most of the lowest gas scale height (by volume) will be
similar to that at the photosphere. Thus dust settling will likely help
prevent thermal waves from suppressing instability.

##  Appendix B Methods for the Model with Realistic Surface

###  B.1 Time Integration

The model equations of motion are listed in §5.2. We evolve the thermal
equation forward in time with the Euler method. The spatial derivative in
equation (2) is, at the $i$th gridpoint, taken to be proportional to
$H/r|_{i+1}-H/r|_{i}$—i.e.,, the first-order upwind scheme. The grid in $r$ is
logarithmic, with typically 250 gridpoints per decade. At each timestep,
$H(r)$ on the spatial grid is obtained from $T(r)$ as follows: after
converting from $T$ to $h$ with the hydrostatic equation, the integral in the
equation for the optical surface (eq. 20) is performed numerically at the
$i$th gridpoint for an assumed value of $H(r_{i})$; and that $H(r_{i})$ is
then adjusted via a root-finder until the integral is unity. For the boundary
conditions, we typically freeze the temperature near the inner and outer
boundary.

###  B.2 Solving for Equilibrium

As stated in §5.1, we obtain the equilibrium disk profile by iteratively
solving the steady state thermal equation (eq. 10) and the photosphere
equation (eq. 20). A naive application of this procedure is typically
numerically unstable. But that instability may be avoided by using the
integral of the thermal equation in place of the thermal equation itself. In
particular, for any assumed $\chi(r)$ profile one may integrate equation (10)
to yield

$\displaystyle{H\over r}=0.02\,\chi_{\rm eff}^{8/7}\,r_{\rm AU}^{2/7}$ (B1)

after using hydrostatic equilibrium and defining $\chi_{\rm eff}$ via ${1\over
r^{2}\chi_{\rm eff}^{8}}\equiv\int_{r}^{\infty}{2d\ln r^{\prime}\over
r^{\prime 2}\chi(r^{\prime})^{8}}$. Therefore the iteration proceeds as
follows: (i) given a profile for $\chi(r)$ (initially taken to be constant),
calculate $H$ via equation (B1), and thence $h=H/\chi$; (ii) use that $h$ in
the equation for the optical surface, and solve for $H(r)$, and thence
$\chi=H/h$; (iii) insert that $\chi$ back into step (i) and repeat until
convergence.

###  B.3 Measuring Linear Growth Rates Numerically

We explain how we measure the exact growth rates shown in Figs. 4 and 7. It is
non-trivial because the growth rate is in general complex-valued. We start
from the linearized thermal equation (eq. 14), which is valid for arbitrary
forms of $F_{\rm irr}[T]$, and take the perturbed temperature to be of the
form ${\delta T\over T}={\rm
Real}(\mbox{\boldmath$\epsilon_{T}$\unboldmath}e^{ik\ln r})$, where in this
appendix we denote complex numbers in bold. As a result, ${\delta F_{\rm
irr}\over F_{\rm irr}}={\rm
Real}(\mbox{\boldmath$\sigma\epsilon_{T}$\unboldmath}e^{ik\ln r})$, where
$\sigma$ is a proportionality constant whose value we determine as follows.
First, we write $\sigma$ and $\epsilon_{T}$ in amplitude-phase form:
$\mbox{\boldmath$\sigma$\unboldmath}=\sigma e^{i\phi_{\sigma}}$ and
$\mbox{\boldmath$\epsilon_{T}$\unboldmath}=\epsilon_{T}e^{i\phi_{T}}$, which
implies that

$\displaystyle{\delta F_{\rm irr}\over F_{\rm
irr}}=\sigma\epsilon_{T}\cos(\phi_{\sigma}+\phi_{T}+k\ln r)\ .$ (B2)

And second, we evaluate the perturbed $F_{\rm irr}$ numerically, via

$\displaystyle{\delta F_{\rm irr}\over F_{\rm irr}}={{F_{\rm
irr}}[T\left(1+{\rm Real}(\mbox{\boldmath$\epsilon_{T}$\unboldmath}e^{ik\ln
r})\right)]\over F_{\rm irr}[T]}-1\,.$ (B3)

To evaluate $\sigma$ at any desired $r$ and $k$, we choose a fixed value of
$\epsilon_{T}\ll 1$ and various values of $\phi_{T}$, and evaluate Equation
(B3) at those values. The resulting function of $\delta F_{\rm irr}/F_{\rm
irr}$ versus $\phi_{T}$ is fit to a sinusoid with the form of equation (B2),
which allows us to extract the values of $\sigma$ and $\phi_{\sigma}$–and
hence $\sigma$.

With the value of $\sigma$ in hand, we substitute into equation (14) the
relation between the complex amplitudes, ${\mbox{\boldmath$\delta F_{\rm
irr}$\unboldmath}\over F_{\rm
irr}}=\mbox{\boldmath$\sigma$\unboldmath}\mbox{\boldmath$\epsilon_{T}$\unboldmath}$.
Equation (14) then yields the solution
$\mbox{\boldmath$\epsilon_{T}$\unboldmath}\propto
e^{\mbox{\boldmath$s$\unboldmath}t}$, where the (complex) growth rate is

$\displaystyle\mbox{\boldmath$s$\unboldmath}={1\over t_{\rm
th}}\left(\mbox{\boldmath$\sigma$\unboldmath}-4\right)\,.$ (B4)

###  B.4 Analytical Growth Rate

We derive the growth rate analytically, under a number of simplifying
assumptions. To begin, we first determine the relationship between the
perturbed temperature and the perturbed heating rate. When the temperature is
perturbed by $\delta T$, hydrostatic equilibrium provides the perturbation in
scale height ($\delta h$), and the equation for the optical surface (eq. 20)
then yields $\delta H$, which in turn leads to the perturbed heating rate.

Under perturbations $\delta h(r)$ and $\delta H(r)$ (but no change in surface
density), the perturbed equation for the optical surface (eq. 20) is Taylor-
expanded to linear order, which yields

$\displaystyle 0=\int_{0}^{r}{dr^{\prime}\over
r^{\prime}}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}\left({h(r^{\prime})\over
r^{\prime}}\right)}e^{-{1\over 2}{H^{2}\over r^{2}}{r^{\prime 2}\over
h(r^{\prime})^{2}}}\times\left\\{{{H^{2}}\over{r^{2}}}\,{{r^{\prime
2}}\over{h(r^{\prime})^{2}}}\,\left[-{{\delta H}\over H}+{{\delta
h(r^{\prime})}\over{h(r^{\prime})}}\right]-{{\delta
h(r^{\prime})}\over{h(r^{\prime})}}\right\\}\,,$ (B5)

where $H$ and $\delta H$ are understood to be functions of $r$. In the limit
that $\chi=H/h\gg 1$, the last term inside the curly brackets is small and we
obtain

$\displaystyle{\delta H\over H}={\rm NUM\over DEN}$ (B6)

where the numerator and the denominator are, respectively,

$\displaystyle{\rm NUM}=\int_{0}^{r}{dr^{\prime}\over
r^{\prime}}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}\left({h(r^{\prime})\over
r^{\prime}}\right)^{3}}{\delta h(r^{\prime})\over h(r^{\prime})}e^{-{1\over
2}{H^{2}\over r^{2}}{r^{\prime 2}\over h(r^{\prime})^{2}}}\,,\hskip
36.135pt{\rm DEN}=\int_{0}^{r}{dr^{\prime}\over
r^{\prime}}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}\left({h(r^{\prime})\over
r^{\prime}}\right)^{3}}e^{-{1\over 2}{H^{2}\over r^{2}}{r^{\prime 2}\over
h(r^{\prime})^{2}}}$ (B7)

One may use a trick to perform the DEN integral: taking the derivative of the
equation for the optical surface (eq. 20), one finds

$\displaystyle 0$ $\displaystyle=$ $\displaystyle{d\over
dr}\int_{0}^{r}dr^{\prime}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}h(r^{\prime})}e^{-{1\over
2}{H(r)^{2}\over r^{2}}{r^{\prime 2}\over h(r^{\prime})^{2}}}$ (B8)
$\displaystyle=$
$\displaystyle{2\tau_{V}(r)\over\sqrt{2\pi}h(r)}e^{-{\chi(r)^{2}\over
2}}+\int_{0}^{r}dr^{\prime}{2\tau_{V}(r^{\prime})\over\sqrt{2\pi}h(r^{\prime})}e^{-{1\over
2}{H(r)^{2}\over r^{2}}{r^{\prime 2}\over h(r^{\prime})^{2}}}\left(-{H(r)\over
r}{d(H/r)\over dr}{r^{\prime 2}\over h(r^{\prime})^{2}}\right)$ (B9)
$\displaystyle=$
$\displaystyle{2\tau_{V}(r)\over\sqrt{2\pi}h(r)}e^{-{\chi(r)^{2}\over
2}}-\left({H\over r}{dH/r\over dr}\right)\times{\rm DEN}$ (B10)

where $\chi(r)=H(r)/h(r)$. The latter equations may be solved for DEN to yield

$\displaystyle{\rm DEN}={1\over{H\over r}{d{H\over r}\over d\ln
r}}{2\tau_{V}\over\sqrt{2\pi}(h/r)}e^{-{\chi(r)^{2}\over
2}}\approx{1\over{\chi^{2}\gamma}}\,{2\tau_{V}\over\sqrt{2\pi}(h/r)^{3}}e^{-{\chi^{2}\over
2}}\,,$ (B11)

after approximating $d\ln{H/r}/d\ln r\approx\gamma$, with $\gamma$ being the
logarithmic slope of the background $h/r$ (eq. 22), as is valid for large
enough $\chi$.

To perform the NUM integral, we write the temperature perturbation to be in
the complex form of eq. (17). We also have ${\delta h\over h}={1\over
2}{\delta T\over T}$ in hydrostatic equilibrium. We may now approximate the
NUM integral as we did for the background approximate equation for the optical
surface (see footnote 3) to obtain

$\displaystyle{\rm NUM}\approx{1\over 2}{1\over ik+\chi^{2}\gamma}{2\tau_{\rm
V}\over\sqrt{2\pi}(h/r)^{3}}e^{-{\chi^{2}\over 2}}\,{{\delta T}\over T}\,\,,$
(B12)

Combining the above two results, we arrive at

$\displaystyle{\delta H\over H}={1\over 2}{\chi^{2}\gamma\over
ik+\chi^{2}\gamma}{\delta T\over T}$ (B13)

For comparison, the simplistic surface model (constant $\chi$) gives $\delta
H/H=\delta h/h={1\over 2}\delta T/T$, which agrees at small $k$.

We may now insert $\delta H/H$ into the heating rate (eq. 2), which yields

$\displaystyle{\delta F_{\rm irr}\over F_{\rm irr}}={{d\over d\ln
r}\left({\delta H\over H}{H\over r}\right)\over{d\over d\ln r}{H\over
r}}=\left({ik\over\gamma}+1\right){\delta H\over H}=\left({\chi^{2}\over
2}{ik+\gamma\over ik+\chi^{2}\gamma}\right)\,{{\delta T}\over T}\ ,$ (B14)

where in the second equality we assume that the spatial variation of Eq. (B13)
is dominated by the sinusoidal dependence, as is reasonable for $k\gtrsim$ a
few. Inserting this into the linearized thermal equation (eq. 14) yields the
complex growth rate

$\displaystyle s\approx{1\over 2t_{\rm
th}}\left[{ik(\chi^{2}-8)-7\chi^{2}\gamma\over ik+\chi^{2}\gamma}\right]\ .$
(B15)

##  Appendix C Numerical Instability in RADMC due to Iterations

As described in §6, a number of papers attempt to find the equilibrium state
of passive disks with an iterative scheme, based on a radiative transfer code
such as RADMC. These often find that the iterations do not converge, but
instead produce large-amplitude propagating waves (e.g., Dullemond & Dominik,
2004b; Min et al., 2009; Siebenmorgen & Heymann, 2012; Ueda et al., 2019).
Here we show that the reason for the non-convergence is an instability closely
connected with the irradiation instability described in this paper. However,
the two instabilities are not identical: the iteration instability is partly
polluted by a numerical (i.e., non physical) instability, and so if the disk
is unstable under iterations, it is not necessarily unstable in reality.

As described in §B.3, the key quantity governing the stability of a disk is
the complex number $\sigma$ that relates heating perturbations to temperature
perturbations, defined via

$\displaystyle{\mbox{\boldmath$\delta{F_{\rm irr}}$\unboldmath}\over F_{\rm
irr}}=\mbox{\boldmath$\sigma$\unboldmath}{\mbox{\boldmath$\delta
T$\unboldmath}\over T}\ .$ (C1)

In the above, the bold $\delta T$ is a complex amplitude, i.e., the real
temperature perturbation is $\delta T={\rm Real}\left(\mbox{\boldmath$\delta
T$\unboldmath}e^{ik\ln r}\right)$, and similarly for $\delta F_{\rm irr}$;
$\sigma$ is a function of both $r$ and $k$, and its value may be determined
numerically (§B.3). Given the value for $\sigma$, the linearized thermal
equation (eq. 14) shows that the disk is unstable if ${\rm
Real}(\mbox{\boldmath$\sigma$\unboldmath})>4$; otherwise, it is stable.

Now, let us compare this behavior with that of the iteration scheme, in which
the temperature at the $k$th iterative step ($T_{k}$) is the equilibrium
solution of the thermal equation, i.e., from eq. (33), $T_{k}=T_{\rm
RADMC}[T_{k-1}]$. Equivalently, writing this in terms of $F_{\rm irr}$ (eq.
32),

$\displaystyle T_{k}=\left({F_{\rm irr}[T_{k-1}]/2\over\sigma_{\rm
SB}}\right)^{1/4}$ (C2)

Setting $T_{k}=T+\delta T_{k}$, where $T$ is the equilibrium temperature, and
linearizing eq. (C2) yields

$\displaystyle{\mbox{\boldmath$\delta T$\unboldmath}_{k}\over T}={1\over
4}{\mbox{\boldmath$\delta F_{\rm irr}$\unboldmath}_{,k-1}\over F_{\rm
irr}}={\mbox{\boldmath$\sigma$\unboldmath}\over 4}{\mbox{\boldmath$\delta
T$\unboldmath}_{k-1}\over T}$ (C3)

The solution to this difference equation is $\mbox{\boldmath$\delta
T$\unboldmath}_{k}={\rm
const}\times\left({\mbox{\boldmath$\sigma$\unboldmath}\over 4}\right)^{k}$.
Therefore the solution is unstable if
$|\mbox{\boldmath$\sigma$\unboldmath}|>4$, which is less stringent than the
true criterion for instability (${\rm
Real}(\mbox{\boldmath$\sigma$\unboldmath})>4$). We note that typically the
real and imaginary parts of $\sigma$ are comparable to each other, and so the
criterion for iterative instabiility is incorrect by an order-unity factor.
For example, in the model with a realistic surface, the marginal $k$ for
stability (plotted in Fig. 4) is $\sim$2 for the iterative instability—rather
than the true value of $\sim 3-4$ as is shown in the figure.

##  Appendix D Gas-Dust Thermal Coupling

In our work, we have assumed that the gas and dust components in the disk
share the same temperature. This requires efficient energy transfer. We
examine the validity of this assumption here.

We first focus on the midplane gas. As gas is inert radiatively, it only
receives or loses energy when colliding with the dust grains. Assuming all
dust particles have size $s$ and bulk density $\rho_{\rm bulk}$, we can
estimate the timescale for changing the gas’s thermal energy as,

$\tau_{\rm gas-dust}\sim{{s\rho_{\rm bulk}}\over{c_{s}\rho_{\rm dust}}}\,.$
(D1)

Comparing this to the dynamical time, and evaluating at the midplane of our
adopted disk (eq. 6), we find

$\Omega\,\tau_{\rm gas-dust}\sim{{s\rho_{\rm bulk}}\over{\Sigma_{\rm
dust}}}\sim 0.005\,\left({s\over{1\rm mm}}\right)\left({{\rho_{\rm
bulk}}\over{1{\rm g/cm^{3}}}}\right)r_{\rm AU}\,.$ (D2)

This process is therefore fast throughout the disk midplane. The strong
coupling between gas and dust may fail at high altitudes. Fortunately, this
should not affect the irradiation instability.