A global two-layer radiative transfer model for axisymmetric, shadowed protoplanetary disks

Okuzumi, Satoshi; Ueda, Takahiro; Turner, Neal J

doi:10.1093/pasj/psac040

Abstract

Understanding the thermal structure of protoplanetary disks is crucial for modeling planet formation and interpreting disk observations. We present a new two-layer radiative transfer model for computing the thermal structure of axisymmetric irradiated disks. Unlike the standard two-layer model, our model accounts for the radial as well as vertical transfer of the starlight reprocessed at the disk surface. The model thus allows us to compute the temperature below “shadowed” surfaces receiving no direct starlight. Thanks to the assumed axisymmetry, the reprocessed starlight flux is given in a one-dimensional integral form that can be computed at a low cost. Furthermore, our model evolves the midplane temperature using a time-dependent energy equation and can therefore treat thermal instabilities. We apply our global two-layer model to disks with a planetary induced gap and confirm that the model reproduces the disks’ temperature profiles obtained from more computationally expensive Monte Carlo radiative transfer calculations to an accuracy of less than 20%. We also apply the model to study the long-term behavior of the thermal wave instability in irradiated disks. Being simple and computationally efficient, the global two-layer model will be suitable for studying the interplay between disks’ thermal evolution and dust evolution.

1 Introduction

Protoplanetary disks’ temperature structure governs many aspects of planet formation. The radial compositional gradient of gas and solids caused by the radial temperature gradient may determine what planets forming at different positions are made of (Öberg et al. 2011). The location of the water snow line, where water ice sublimates and condenses, is particularly important in this context as it may constrain where rocky planets like the Earth form. Ice’s sublimation, condensation, and sintering around the snow line can also cause a change in solid particles’ stickiness (e.g., Chokshi et al. 1993; Dominik & Tielens 1997; Wada et al. 2009; Sirono & Ueno 2017)¹ and a local pileup of silicates and ice (Stevenson & Lunine 1988; Cuzzi & Zahnle 2004; Saito & Sirono 2011; Sirono 2011; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Hyodo et al. 2019). These processes may not only affect planet formation directly but may also produce some observable features in disk thermal emission (Banzatti et al. 2015; Zhang et al. 2015; Okuzumi et al. 2016; Pinilla et al. 2017). Knowledge of disks’ thermal structure is also necessary for interpreting scattered light’s radial falloff and planet-launched spiral waves (e.g., Isella & Turner 2018), inferring turbulence from super-thermal linewidths (e.g., Flaherty et al. 2018, 2020), and guiding planets’ orbital migration (e.g., Paardekooper et al. 2010, 2011; Bitsch et al. 2013).

Despite its importance, our understanding of disks’ thermal structure is still limited. This is particularly true for the thermal structure deep inside the disks, where planet formation mainly occurs. The disk thermal structure well above the midplane has been well studied with observations of spectral energy distribution (e.g., Kenyon & Hartmann 1987; Chiang & Goldreich 1997; D’Alessio et al. 1998, 1999; Chiang et al. 2001; Sierra & Lizano 2020) and optically thick molecular emission lines (e.g., Akiyama et al. 2011; Rosenfeld et al. 2013; Weaver et al. 2018; Calahan et al. 2021). The temperature structure closer to the midplane can also be constrained from intensity maps of marginally optically thin emission lines (Zhang et al. 2017) and channel maps of optically thick lines (Dullemond et al. 2020), but there are only a few cases for which these approaches have been applied. The midplane temperature profile can also be inferred from multiwavelength imaging of dust continuum emission (Kim et al. 2019; Carrasco-González et al. 2019; Macías et al. 2021), but the inferred temperature depends on the assumed scattering properties of the opacity-dominating dust grains (see Carrasco-González et al. 2019).

The disk thermal structure is determined by stellar radiation, internal heating associated with disk accretion, and radiative cooling by gas and dust. In the classical framework of the viscous accretion model (Lynden-Bell & Pringle 1974), accretion heating is the dominant heating mechanism in the inner few au midplane region, thus controlling the location of the snow line (e.g., Davis 2005; Garaud & Lin 2007; Oka et al. 2011; Bitsch et al. 2015). However, recent studies of the magnetohydrodynamics of weakly ionized protoplanetary disks have shown that accretion heating favors the disk surface because the midplane region is too poorly ionized to sustain a large electric current (Hirose & Turner 2011; Mori et al. 2019). These models predict that magnetohydrodynamic accretion heating can raise the midplane temperature only when the disk opacity is high enough (Béthune & Latter 2020; Mori et al. 2021). This suggests that stellar irradiation rather than internal heating may determine the location of the water snow line.

The study of the thermal structure of passively irradiated protoplanetary disks has a long history (e.g., Kusaka et al. 1970; Calvet et al. 1991; Chiang & Goldreich 1997; D’Alessio et al. 1998; Dullemond et al. 2001, 2002), but its complete understanding has yet to be established because of two complications. The first complication is that the disk thermal structure entirely depends on the global distribution of small dust grains that absorb and reprocess stellar radiation. Because starlight grazes the disk surface at a small angle (typically ∼0.01–0.1 rad), even a small perturbation on the surface can yield a shadowed region that does not directly receive stellar radiation (Dullemond et al. 2001; Dullemond & Dominik 2004a, 2004b). Shadowed regions can have significantly low temperatures and hence can have gas and dust compositions that differ greatly from their surroundings (Ohno & Ueda 2021). However, to accurately compute the temperature of shadowed regions, one must account for the radial radiative transfer of reprocessed starlight (e.g., Jang-Condell & Sasselov 2003; Dullemond & Dominik 2004a; Turner et al. 2012; Jang-Condell & Turner 2012; Ueda et al. 2017). The growth, settling, and radial migration of small dust grains should also be taken into account as these processes may change the distribution of the shadows.

The second complication concerns the instabilities of disks’ thermal structure. Dullemond (2000) and Watanabe and Lin (2008) showed that optically thick disks irradiated by the central star are unstable to self-shadowing in the limits of short and long thermal relaxation timescales, respectively (for more recent studies, see Siebenmorgen & Heymann 2012; Ueda et al. 2021; Wu & Lithwick 2021; Pavlyuchenkov et al. 2022). This instability, which we call the thermal wave instability after Watanabe and Lin (2008)², is intrinsically related to the starlight’s small grazing angle mentioned above: even a small hill on the irradiation surface can cast a shadow. The hill’s starlit inner side warms and expands, while its shadowed outer side cools and contracts, so that the hill propagates towards the star (see figure 1 of Ueda et al. 2021 and Wu & Lithwick 2021 for a cartoon describing the mechanism of the instability). The surface waves generated by the instability produce the interior temperature’s fluctuations that propagate inward on a timescale comparable to the local thermal relaxation timescale at the midplane. The instability is potentially relevant to planet formation because it can cause snow lines to oscillate radially (Ueda et al. 2021) and may even produce pressure maxima collecting pebble-sized particles (Watanabe & Lin 2008). More recently, Owen (2020) has identified a distinct type of thermal instability that results from an abrupt change in the opacity across a snow line. However, fully understanding these thermal instabilities requires a radiative transfer model that can treat shadowing and does not assume radiative equilibrium. Lacking a simple model that fulfills these requirements, the roles of the thermal instabilities in planet formation have not been elucidated so far.

Fig. 1.

Schematic showing the transfer of stellar radiation in (a) a uniformly flared disk and in (b) a flared disk with an annular gap, i.e., a dip in the radial surface density profile. The disk surface here is the location where vertical optical depth at visible wavelengths reaches ∼H/r, where H is the disk’s scale height. For flared disks, this surface well approximates the true irradiation surface where the optical depth along starlight rays reaches unity (section 2). Flared disk surfaces receive direct (visible) starlight and reprocess it into infrared (IR) radiation. The conventional, local two-layer model computes the disk’s interior temperature using the vertical transfer of the reprocessed starlight (a). The local two-layer model does not apply to deep gaps that receive no direct starlight. To obtain the temperature in such a shadowed gap, one must consider the radial transfer of reprocessed starlight (b).

Open in new tab Download slide

In this study, we present a simple radiative transfer model that can treat shadowed protoplanetary disks. Our model is a generalization of the well-known two-layer model (Kusaka et al. 1970; Chiang & Goldreich 1997; Chiang et al. 2001), which computes the interior temperature of an optically thick disk by considering the starlight reprocessed on the disk surface as the heating source. Being extremely simple and moderately accurate (see Dullemond & Dominik 2004a), the two-layer model has been widely used to model the radial temperature distribution of irradiated protoplanetary disks. However, the conventional two-layer model only considers the vertical transfer of the reprocessed starlight and hence is unable to compute temperatures just below shadowed surfaces. Our new model, which we call the global two-layer model, resolves this issue by accounting for the radial transfer of the reprocessed starlight. Our model is limited to axisymmetric disks, but is much simpler than exising radiative transfer models that can treat shadows, including those based on the Monte Carlo approach (Dullemond et al. 2012; Turner et al. 2012) and those considering three-dimensional perturbations on the disk surface (Jang-Condell & Sasselov 2003, 2004; Jang-Condell 2008, 2009; Jang-Condell & Turner 2012). Therefore, our model will be suitable for coupled simulations of dust evolution and disk thermal evolution. Furthermore, our model does not rely on radiative equilibrium and hence can treat disks’ thermal instabilities. The aims of this paper are to formulate the global two-layer model and demonstrate its applicability to shadowed disks.

The structure of this paper is as follows. Section 2 reviews the standard local two-layer model relying on the plane-parallel approximation and highlights its limitation. Section 3 describes our new two-layer model, and section 4 tests the model by comparing it with Monte Carlo radiative transfer calculations for gapped disks. In section 5, we use the model to simulate long-term evolution of the thermal wave instability. Section 6 presents a summary and future directions.

2 Local two-layer model

We begin by reviewing the conventional local two-layer radiative transfer model for passively irradiated disks (Chiang & Goldreich 1997; Chiang et al. 2001; see also Krügel 2008; Armitage 2010 for textbook descriptions), with a particular emphasis on its key assumptions and limitations. The aim here is to explain why the local two-layer model is inapplicable to disks with shadows.

The local two-layer model assumes that (1) the disk is locally plane-parallel, i.e., the disk structure changes radially on a scale much longer than the disk vertical thickness, and that (2) the disk is flared, so that every portion of the disk surface receives direct starlight (see figure 1 for a schematic). Here, the disk surface refers to where the vertical optical depth at optical wavelengths is ∼H/r, where H is the disk’s scale height. For flared disks, this surface well approximates the true irradiation surface where the optical depth to the grazing starlight reaches unity because the radial optical depth is ∼r/H times larger than the vertical depth (in other words, the starlight grazing angle for flared disks is comparable to H/r; see, e.g., Chiang & Goldreich 1997).

With the plane-parallel and flared assumptions, the vertical downward flux of the reprocessed (infrared) radiation from the irradiation surface toward the disk interior is given by

$$\begin{equation} F_{\rm rep, \downarrow } = f_\downarrow \mu _* F_*, \end{equation}$$

(1)

where F_* is the magnitude of the direct starlight flux, μ_* is the sine of the grazing angle between the starlight and disk surface, and f_↓ a dimensionless coefficient that depends on the single-scattering albedo of the grains receiving the starlight. In equation (1), the product μ_*F_* represents the flux of the stellar radiation incident on the irradiation surface, while the prefactor f_↓ represents the fraction of the starlight flux reprocessed downward. A purely absorbing atmosphere gives f_↓ = 1/2, meaning that it reprocesses half of the incident stellar radiation energy into downward thermal radiation (Chiang & Goldreich 1997). For a plane-parallel atmosphere of nonzero albedo, one has (see appendix 1 for a derivation)

$$\begin{equation} f_\downarrow = \frac{\epsilon _*}{2} + \left(\sqrt{3\epsilon _*}+\frac{3}{2q}\right)\frac{1-\epsilon _*}{3+2\sqrt{3\epsilon _*}}, \end{equation}$$

(2)

where

$$\begin{equation} \epsilon _* = \frac{\kappa _{*} }{\chi _{*} } \end{equation}$$

(3)

is the ratio between the Planck-mean absorption and extinction opacities κ_* and χ_* evaluated at the stellar surface temperature T_*, and

$$\begin{equation} q = \frac{\chi _*}{\chi _{\rm R}} \end{equation}$$

(4)

is the ratio between χ_* and the Rosseland mean extinction opacity for the disk’s own thermal radiation, χ_R [see equations (A3), (A4), and (A8) for a definition of the mean opacities]. The ratio ε_* is related to the grains’ single-scattering albedo ω_* for the starlight as ω_* = 1 − ε_*. The fraction f_↓ decreases monotonically with decreasing ε_*, reflecting the fact that multiple scattering by the gas and dust particles enhances the backscattering of the starlight. Figure 2 illustrates f_↓ as a function of ε_* for the particular case of q = 5.7 considered in section 4.

Fig. 2.

Fraction f_↓ of the stellar radiation reprocessed downward as a function of ε_* for the case q = χ_*/χ_R = 5.7 considered in section 4. Here, ε_* is unity minus the grains’ starlight albedo. A value of ε_* < 1 yields f_↓ < 1/2; for instance, ε_* = 0.12 gives f_↓ = 0.24.

Open in new tab Download slide

Assuming that the disk interior is optically thick to infrared radiation and is also in radiative equilibrium, the downward infrared flux given by equation (1) balances the upward thermal flux from the interior, F_i = σ_SBT⁴, where T is the interior temperature and σ_SB is the Stefan–Boltzmann constant. This balance yields

$$\begin{equation} T(r) = \biggl ({ \displaystyle \frac{f_\downarrow \mu _* F_*}{\sigma _{\rm SB}}} \biggr ) ^{1/4} . \end{equation}$$

(5)

For a central star of radius R_* and surface temperature T_*, F_* can be written as

$$\begin{equation} F_* = \frac{\sigma _{\rm SB} R_*^2 T_*^4}{r^2+z_{\rm s}(r)^2} \approx \frac{\sigma _{\rm SB} R_*^2 T_*^4}{r^2}, \end{equation}$$

(6)

where r is the cylindrical distance from the central star and z_s(r) is the height of the irradiation surface at r. The final expression assumes z_s(r)²/r² ≪ 1, which holds in typical protoplanetary disks. Substituting equation (6) into equation (5) yields

$$\begin{equation} T(r) = (f_\downarrow \mu _*)^{1/4} \biggl ({ \displaystyle \frac{R_*}{r}} \biggr ) ^{1/2} T_*. \end{equation}$$

(7)

For purely absorbing atmospheres (f_↓ = 1/2), equation (7) reduces to the well-known expression for the disk interior temperature in the conventional two-layer model neglecting scattering [e.g., equation (4) of Dullemond et al. 2001.

It is important to note that every quantity in the right-hand side of equation (5), or equivalently of equation (7), is evaluated locally. In other words, the conventional two-layer model only accounts for the vertical transfer of reprocessed starlight as schematically shown in figure 1a.

However, such a model is incompatible with disks with a shadow. To illustrate this, we consider a disk with an annular gap, i.e., a local dip in the radial surface density profile (figure 1b). For the moment, we continue to refer to the location where the vertical visible optical depth is ∼H/r as the disk surface. As shown in figure 1b, this surface has a valley in the gap. We now assume that the valley is so deep that it falls into the shadow cast by the region inward of the gap. The problem is that such a valley receives no direct starlight and therefore provides no reprocessed radiation toward the midplane. Hence, one cannot use equations (5) or (7) to compute the midplane temperature in the shadowed gap. This example clearly indicates that the radial transfer of reprocessed radiation must be taken into account to compute the temperature in shadowed regions (see figure 1b).

3 Global two-layer model

In this section, we present our global two-layer model that takes into account the radial as well as vertical transfer of reprocessed starlight. As in the standard two-layer model, we consider the irradiation surface reprocessing direct (visible) starlight into thermal (infrared) radiation and the disk interior heated by the reprocessed starlight (subsections 3.1 and 3.2). We relax the two fundamental assumptions of the standard two-layer model by splitting the irradiation surface into concentric rings and calculate the two-dimensional transfer of the downward emission from individual rings (subsection 3.3). Our current model assumes vertical hydrostatic equilibrium (subsection 3.4) but does not assume radiative equilibrium (subsection 3.5). We also present numerical implementation of the model (section 3.6).

Our global two-layer model is largely inspired by the radiative transfer model of Jang-Condell and Sasselov (2003, 2004) and Jang-Condell (2008), which splits the irradiation surface into small elements and computes three-dimensional radiative transfer of the reprocessed starlight from each element. Compared to this previous model, our model is limited to axisymmetric disks. However, this assumed symmetry allows us to greatly simplify the expression of the reprocessed starlight flux (for details, see subsection 3.3).

In the following, we employ the cylindrical coordinate system (r, φ, z) centered on the central star. Because we assume axisymmetry, every quantity is independent of φ.

3.1 The irradiation surface

To properly treat shadows, we do not use an approximate irradiation surface defined in terms of the vertical optical depth as used in figure 1. We employ the exact irradiation surface defined as where the optical depth to the direct starlight is unity (D’Alessio et al. 2001). At position (r, z) in a disk, the Planck-mean optical depth to the direct starlight is given by

$$\begin{equation} \tau _{*} = \int _{\rm starlight} \chi _{*} \rho \, ds, \end{equation}$$

(8)

where χ_* is the Planck-mean extinction opacity per gas mass for the starlight, ρ is the gas mass density, and s denotes the path length of the stellar ray from the star to the position considered. Thus, the height of the irradiation surface, z_s(r), at arbitrary cylindrical radius r is defined by the relation τ_*[r, z_s(r)] = 1. The disk’s vertical structure is specified in subsection 3.4.

At large radial distances where the central star can be seen as a point source, stellar rays travel along lines of approximately constant z/r. For such regions, it is useful to map the irradiation surface in the r–z/r plane. Figure 3 schematically shows the irradiation surface of the gapped disk shown in figure 1b. As already mentioned in section 2, the irradiation surface outside the shadowed gap approximately matches the surface at which the vertical visual optical depth is ∼H/r. Inside the shadowed gap, the irradiation surface as seen in the r–z/r plane is represented by a horizontal line. Below this line, the material inward of the gap blocks the direct starlight (see the lower panel of figure 3).

Fig. 3.

Upper panel: Schematic showing a disk’s irradiation surface (dotted curve; see subsection 3.1) and thermal (infrared) photosphere (thin solid curve; see subsection 3.2). The quantities α_*, z_s, and z_IR stand for the starlight grazing angle, irradiation surface’s height, and thermal photosphere’s height, respectively. Lower panel: Irradiation surface as seen in the r–z/r plane. In this plane, stellar rays of different propagation angles are represented by horizontal lines. A shadowed region refers to where α_* vanishes. As depicted in figure 1b, this region receives no direct stellar radiation but is heated by the reprocessed starlight that propages both radially and vertically.

Open in new tab Download slide

For given z_s(r), one can calculate μ_* = sin α_* as (Kusaka et al. 1970; Ruden & Pollack 1991)

$$\begin{equation} \mu _* = \sin ^{-1} \biggl ({ \displaystyle \frac{4R_*}{3\pi r}} \biggr ) + \tan ^{-1}\left(\frac{z_{\rm s}}{r}\frac{d\ln z_{\rm s}}{d\ln r}\right) - \tan ^{-1}\left(\frac{z_{\rm s}}{r}\right) . \end{equation}$$

(9)

In the right-hand side of equation (9), the first term accounts for the finite stellar size, while the sum of the second and third terms is related to the grazing angle of the ray from the central point source. At large radial distances where the first term is negligible, μ_* vanishes in shadowed regions with radially constant z_s/r (see the lower panel of figure 3).

3.2 The thermal photosphere

We define the thermal photosphere as the surface on which the optical depth to reprocessed starlight from the irradiation surface reaches unity (see figure 3). Strictly speaking, the location of the thermal photosphere depends on the incident angle of the reprocessed radiation. To avoid this complexity, we approximate the thermal photosphere by the surface where the vertical optical depth for downward infrared radiation exceeds 1/2. This choice is based on the fact that the angle-averaged optical depth of an optically thin disk to isotropic radiation is twice the vertical optical depth (Nakamoto & Nakagawa 1994). Writing the Planck-mean extinction opacity for infrared radiation as χ_P, the corresponding vertical optical depth is given by

$$\begin{equation} \tau _{\rm P}(r,z) = \int ^\infty _z \chi _{\rm P}(r,z^{\prime })\rho (r,z^{\prime })\, dz^{\prime }. \end{equation}$$

(10)

Using this, the height of the thermal photosphere, z_IR(r), at radial position r can be approximately defined by the relation τ_P[r, z_IR(r)] = 1/2. The region well below the thermal photosphere is also optically thick to its own thermal emission.

For disks with τ_P < 1/2 at the midplane, we set z_IR = 0 for convenience. However, the midplane of such an optically thin disk is not a real thermal photosphere, absorbing only a small fraction of the incoming infrared radiation. The low absorptivity and emissivity of optically thin disks are taken into account in subsection 3.5.

3.3 Flux of the downward reprocessed starlight

We split the irradiation surface into thin concentric ring elements of radius r′, radial extent dr′, and height z_s(r′) (figure 4), Each ring emits infrared radiation with a luminosity proportional to the received stellar flux. To derive an analytic expression for the downward radiation flux from each ring element, we further divide the rings azimuthally into rectangular segments of azimuthal width r′dφ′. The area of each segment is |$dS = \sqrt{dr^{\prime 2} + dz_{\rm s}^2} \, r^{\prime } d\phi ^{\prime } = \sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } dr^{\prime } d\phi ^{\prime }$|⁠. The luminosity of the downward radiation from each segment is

$$\begin{eqnarray} dL_{\rm rep, \downarrow } &=& f_\downarrow (r^{\prime }) \mu _* F_* dS \nonumber \\ &=& f_\downarrow \mu _* F_* \sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } d r^{\prime } d\phi ^{\prime }. \end{eqnarray}$$

(11)

The radiation from each surface segment strikes the thermal photosphere. For simplicity, we here assume dz_IR/dr ≪ 1 and approximate the thermal photosphere as a plane locally parallel to the midplane. Because of this approximation, our model neglects the shadowing of the reprocessed radiation field on the thermal photosphere. For typical protoplanetary disks with z_IR ≪ r, the approximation made here can be justified unless z_IR varies on a radial length scale ≪r. Because we assume axisymmetry, the azimuthal coordinate of an arbitrary segment of the thermal photosphere can be taken to be zero without loss of generality. The angle between the reprocessed starlight and the thermal photosphere segment can then be written by

$$\begin{equation} \alpha _{\rm rep} = \sin ^{-1} \biggl ({ \displaystyle \frac{h}{\Delta }} \biggr ) , \end{equation}$$

(12)

where

$$\begin{equation} h(r,r^{\prime }) \equiv z_{\rm s}(r^{\prime }) - z_{\rm IR}(r) \end{equation}$$

(13)

is the height of the irradiation surface segment relative to the photosphere segment and

$$\begin{equation} \Delta \equiv \sqrt{(r-r^{\prime } \cos \phi ^{\prime } )^2 + r^{\prime 2} \sin ^2\phi ^{\prime } + h^2} \end{equation}$$

(14)

is the distance between the two segments.

$Schematic showing how the global two-layer model computes the flux of the reprocessed starlight from the irradiation surface to the optically thick interior. The irradiation surface is decomposed into thin ring elements of radius r′ lying at height z = zs(r′) above the midplane. Each ring element is further divided into segments of azimuthal width r′dφ′. The disk interior that is optically thick to infrared thermal radiation receives the reprocessed starlight on its surface (thermal photosphere) lying at height z = zIR(r). The flux of the reprocessed starlight from each ring element to position [r, zIR(r)] is calculated by integrating the contributions $d{\boldsymbol F}_{\rm rep,\downarrow }$ from the constituting segments.$

Fig. 4.

Schematic showing how the global two-layer model computes the flux of the reprocessed starlight from the irradiation surface to the optically thick interior. The irradiation surface is decomposed into thin ring elements of radius r′ lying at height z = z_s(r′) above the midplane. Each ring element is further divided into segments of azimuthal width r′dφ′. The disk interior that is optically thick to infrared thermal radiation receives the reprocessed starlight on its surface (thermal photosphere) lying at height z = z_IR(r). The flux of the reprocessed starlight from each ring element to position [r, z_IR(r)] is calculated by integrating the contributions |$d{\boldsymbol F}_{\rm rep,\downarrow }$| from the constituting segments.

Open in new tab Download slide

Assuming that the downward reprocessed starlight is isotropic, and approximating the irradiation surface segment as a point source, the magnitude of the downward reprocessed starlight flux |$d{\boldsymbol F}_{\rm rep, \downarrow }$| (see figure 4) can be written as

$$\begin{equation} |d{\boldsymbol F}_{\rm rep, \downarrow }| = \frac{dL_{\rm rep, \downarrow }}{2\pi \Delta ^2}, \end{equation}$$

(15)

where the factor 2π comes from the solid angle of the lower hemisphere. Because the reprocessed starlight strikes at angle α_rep, its flux incident on the thermal photosphere is |$\sin \alpha _{\rm rep}|d{\boldsymbol F}_{\rm rep, \downarrow }|$|⁠. Integrating this over φ′, we obtain

$$\begin{eqnarray} && \int _{\phi ^{\prime }=0}^{\phi ^{\prime }=2\pi } \sin \alpha _{\rm rep}|d{\boldsymbol F}_{\rm rep, \downarrow }| \nonumber \\ &&\quad= \frac{f_\downarrow \mu _* F_* h\sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } dr^{\prime } }{2\pi } \int _0^{2\pi } \frac{d\phi ^{\prime } }{\Delta ^3} \nonumber \\ &&\quad= \frac{2f_\downarrow \mu _* F_* h \sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } d r^{\prime }\, E(k)}{\pi [ (r^{\prime }-r)^2 + h^2] \sqrt{(r^{\prime }+r)^2 + h^2} }, \end{eqnarray}$$

(16)

where |$E(k) = \int _0^{\pi /2} (1-k^2 \sin ^2 x)^{1/2} dx$| is the complete elliptic integral of the second kind with |$k =2\sqrt{rr^{\prime }/[(r^{\prime }+r)^2+h^2]}$|⁠. Finally, by integrating equation (16) over r′, we obtain the net downward vertical flux F_{rep, ↓}(r) of the reprocessed starlight from the entire irradiation surface to the disk interior at radius r,

$$\begin{eqnarray} F_{\rm rep, \downarrow }(r) &=& \int _{\phi ^{\prime }, r^{\prime }} \sin \alpha _{\rm rep}|d{\boldsymbol F}_{\rm rep, \downarrow }| \nonumber \\ &=& \frac{2}{\pi } \int \frac{f_\downarrow \mu _* F_* h\sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } E(k)}{\left[ (r^{\prime }-r)^2 +h^2\right] \sqrt{(r^{\prime }+r)^2 +h^2} }dr^{\prime }. \nonumber \\ \end{eqnarray}$$

(17)

Here, the integration is performed over regions of α_rep > 0 (h > 0)³. Note that f_↓, μ_*, F_*, and h generally depend on r′ and hence should be inside the integral.

Equation (17) is a generalization of the reprocessed starlight flux in the local two-layer model, equation (1). To see this, we rewrite equation (17) as

$$\begin{eqnarray} F_{\rm rep, \downarrow }(r) = \int f_\downarrow (r^{\prime }) \mu _*(r^{\prime })F_*(r^{\prime }) W(r,r^{\prime })dr^{\prime }, \end{eqnarray}$$

(18)

where the function W(r, r′) is

$$\begin{equation} W(r,r^{\prime }) = \frac{2h\sqrt{1 + (dz_{\rm s}/dr^{\prime })^2}\, r^{\prime } E(k)}{\pi \left[ (r^{\prime }-r)^2 +h^2\right] \sqrt{(r^{\prime }+r)^2 +h^2} } \end{equation}$$

(19)

for h > 0 and W(r, r′) = 0 otherwise. Equation (18) can be viewed as the radial average of equation (1) weighted by W(r, r′). As shown in appendix 2, equation (18) recovers equation (1) in the limit of small h.

The weighting function W represents the radially nonlocal nature of the reprocessed radiation. As explained below, this function is peaked around r′ = r and has a width of ∼h. Because the elliptic integral E is bounded in the narrow range of 1 ≤ E ≤ π/2 and because |$\sqrt{1 + (dz_{\rm s}/dr^{\prime })^2} \sim 1$|⁠, the profile of W is essentially determined by the factor |$hr^{\prime }/ \lbrace [ (r^{\prime }-r)^2 +h^2] \sqrt{(r^{\prime }+r)^2 +h^2} \rbrace$|⁠. Since h ≪ r, the profile around the peak is approximately Lorentzian ∝ 1/((r − r′)² + h²) with the half width at half maximum of h. Figure 5 illustrates the functional form of W for z_s(r′) = 0.15r′ and z_IR(r) = 0.05r, i.e., h = 0.15r′ − 0.05r. In this particular example, W vanishes at r′/r < 1/3, at which h becomes negative.

Fig. 5.

Weighting function W(r, r′) in the radial integration in equation (18) as a function of r′/r for z_s = 0.15r′ and z_IR = 0.05r, normalized such that W = 1 at r′ = r.

Open in new tab Download slide

3.4 Disk vertical structure

So far, we have not specified the disk’s vertical structure. In the present study, we avoid detailed modeling of the vertical structure by approximating the disk as vertically isothermal and hydrostatic. These two approximations yield the vertical distribution of ρ of the simple analytic form

$$\begin{equation} \rho (r,z) = \frac{\Sigma (r)}{\sqrt{2\pi } H(r)}\exp \left[-\frac{z^2}{2H(r)^2}\right], \end{equation}$$

(20)

where Σ is the gas surface density and H = c_s/Ω_K is the gas scale height, with c_s and Ω_K being the isothermal sound speed and Keplerian frequency, respectively. The isothermal sound speed is related to the disk interior temperature as |$c_s = \sqrt{k_{\rm B}T/m_{\rm g}}$|⁠, where k_B is the Boltzmann constant and m_g is the mean mass of the gas molecules.

We comment on the validity and limitations of the vertical isothermal and hydrostatic approximations. For passively irradiated disks in radiative equilibrium, the vertical isothermal approximation is valid well below the irradiation surface. However, this approximation breaks down across the irradiation surface, above which the temperature is higher because of direct stellar irradiation (Calvet et al. 1991; Chiang & Goldreich 1997; D’Alessio et al. 1998). A more reasonable approximation would be to assign different temperatures to the interior (z < z_s) and warmer surface region (z > z_s) as described in appendix B of Watanabe and Lin (2008). However, such two-temperature modeling would result in a vertical density profile dependent on z_s, requiring iterative determination of z_s and ρ(z). We defer this complication to future work. Regarding the equilibrium disk structure, our vertical isothermal approximation should suffice to determine z_s because the warmer but optically thin gas well above the irradiation surface should have little effect on the location of the irradiation surface (see Watanabe & Lin 2008). Nevertheless, the vertical variation of the temperature, and more importantly of the thermal relaxation time, could have a more critical impact on the time evolution of the disk temperature. We discuss this point in more detail in subsection 3.5.

The vertical hydrostatic approximation is valid as long as the disk’s thermal relaxation time t_th [see equation (27) for a definition] is much longer than the orbital timescale t_K = 2π/Ω_K. In protoplanetary disks, the condition t_th ≫ t_K is typically fulfilled at |$r \ll 100\, \rm au$| (Dullemond 2000; Watanabe & Lin 2008; Wu & Lithwick 2021; see also subsection 5.2) but can break down farther out. As discussed by Wu and Lithwick (2021), the disk’s hydrodynamic response may influence its thermal and dynamical evolution. In this study, we concentrate on the inner disk region where the vertical hydrostatic assumption is applicable.

The opacities used in our model depend on the dust-to-gas mass ratio and size distribution of the opacity-dominating dust grains. One can account for dust settling by using vertically varying opacities.

3.5 Disk temperature evolution

Because passively irradiated disks can be thermally unstable, we treat the disk interior temperature as intrinsically time-dependent. For simplicity, we neglect radial heat advection by accreting gas (see, e.g., Cannizzo 1993), radial diffusion of the disk’s own thermal radiation (Latter & Balbus 2012; Owen & Armitage 2014), and accretion heating. We plan to include these effects in future work.

Neglecting the above-mentioned effects, the time evolution of the interior temperature T obeys the following energy equation:

$$\begin{equation} \frac{\gamma +1}{2(\gamma -1)}\frac{k_{\rm B}\Sigma }{m_{\rm g}}\frac{\partial T}{\partial t} = 2{\cal C}[F_{\rm rep, \downarrow } + \sigma _{\rm SB}(T_{\rm ex}^4 - T^4)], \end{equation}$$

(21)

where γ is the adiabatic index (taken to be 1.4 throughout this paper), T_ex is the radiation temperature of the parent molecular cloud (Ueda et al. 2021), and |${\cal C}$| is a dimensionless factor correcting for the effects of the disk’s infrared optical thickness and albedo.

Equation (21) is the vertically integrated equation of total energy conservation (Watanabe et al. 1990; Watanabe & Lin 2008). Its left-hand side stands for the rate of change in the total energy per unit disk area. The temperature on the left-hand side originally stands for the density-weighted vertical average of the temperature (Watanabe et al. 1990). Following Watanabe and Lin (2008), we have applied the vertical isothermal approximation and represented the vertically averaged temperature with the single interior temperature. In reality, the rate of change in temperature in optically thick disks depends on the depth from the disk surface, with shallower and optically thinner regions having shorter cooling timescales. Recently, Pavlyuchenkov, Maksimova, and Akimkin (2022) have pointed out that the vertical variation of the thermal relaxation timescales could weaken the thermal wave instability around the midplane of optically thick disks. This effect is not included in our current modeling.

The right-hand side of equation (21) accounts for radiative heating and cooling on both sides (z < 0 and z > 0) of the disk. We take the correction factor |${\cal C}$| to be

$$\begin{equation} {\cal C} = \frac{4 \sqrt{\epsilon }\tanh (\tau _{\rm eff,mid})}{\sqrt{3} +2\sqrt{\epsilon } \tanh (\tau _{\rm eff,mid})} \frac{1}{1+3\tau _{\rm R,mid}/4}, \end{equation}$$

(22)

where

$$\begin{equation} \epsilon = \frac{\kappa _{\rm P}}{\chi _{\rm R}} \end{equation}$$

(23)

is the ratio between the Planck-mean absorption opacity κ_P [see equation (A7) for a definition] and Rosseland extinction opacity χ_R,

$$\begin{equation} \tau _{\rm eff,\rm mid} =\int _0^\infty \sqrt{3\epsilon }\, \chi _{\rm R}\rho \, dz^{\prime } =\int _0^\infty \sqrt{3\chi _{\rm R}\kappa _{\rm P}}\, \rho \, dz^{\prime } \end{equation}$$

(24)

is the effective absorption optical depth to the midplane accounting for multiple scattering (Rybicki & Lightman 1979), and

$$\begin{equation} \tau _{\rm R,mid} = \int _0^\infty \chi _{\rm R}\rho \, dz^{\prime } \end{equation}$$

(25)

is the Rosseland-mean vertical optical depth to the midplane. All the mean opacities appearing here are for the disk thermal radiation. In equation (22), the factor involving τ_{eff, mid} corrects for the infrared emissivity and absorptivity of both optically thick and thin disks, derived by approximating a disk with a uniform slab [see equation (A21) in appendix 4]. The factor (1 + 3τ_{R, mid}/4)⁻¹ corrects for the vertical radiative diffusion flux in optically thick (τ_{R, mid} ≫ 1) disks being inversely proportional to τ_{R, mid} (Wu & Lithwick 2021).

3.6 Numerical implementation

We consider a radial computational domain spanning r = r_min to r = r_max and discretize it into logarithmically spaced cells. Each radial cell has inner and outer boundaries at r_{j − 1/2} = r_min(r_max/r_min)^{(j − 1)/J} and r_{j + 1/2} = r_min(r_max/r_min)^j/J, respectively, and a logarithmic center at |$r_j = \sqrt{r_{j-1/2}r_{j+1/2}} = r_{\rm min} (r_{\rm max}/r_{\rm min})^{(j-1/2)/J}$|⁠, where j = 1, 2, …, J labels the cells and J is the number of the cells. The cell width should be sufficiently smaller than the width ∼h of the weighting function W.

The irradiation surface is found with a ray-tracing approach. We use rays emanating from the coordinate origin (r, z) = (0, 0) at angle θ with respect to the r-axis. All simulations presented in this paper adopt 180 linearly spaced θ grids spanning θ = 0 to θ = θ_max = π/12, with a resolution of dθ = θ_max/180 = 1.45 × 10⁻³. We have also tried simulations with dθ = θ_max/360 and confirmed that the higher angular resolution gives no appreciable change in the simulation results.

When searching for the irradiation surface, one must assume the starlight optical depth to the inner computational boundary. The starlight optical depth τ_* [equation (8)] to arbitrary radial position r along a ray with angle θ can be written as

$$\begin{equation} \tau _* = \tau _{*,\rm in} + \frac{1}{\cos \theta }\int _{r_{\rm min}}^r \chi _* \rho (r^{\prime },r^{\prime } \tan \theta ) dr^{\prime }, \end{equation}$$

(26)

where |$\tau _{*,\rm in}$| denotes the optical depth to the inner computational boundary. The second term in the right-hand side of equation (26) uses z = rtan θ and ds = dr/cos θ. The problem here is that |$\tau _{*,\rm in}$| is intrinsically unknown because it is determined by the disk structure outside the computational domain. However, as demonstrated in appendix 3, an unreasonable choice of |$\tau _{*,\rm in}$| can cause unwanted artifacts in the resulting temperature distribution near the inner boundary. Our prescription for |$\tau _{*,\rm in}$| is presented in appendix 3.

The energy equation (21) is solved with a first-order forward differencing scheme. The time step must be smaller than the thermal relaxation timescale for the vertically averaged temperature (Watanabe & Lin 2008):

$$\begin{equation} t_{\rm th} = \frac{\gamma +1}{2(\gamma -1)}\frac{k_{\rm B}\Sigma }{m_{\rm g} \sigma _{\rm SB}T^3{\cal C}}. \end{equation}$$

(27)

Note that t_th depends on the correction factor |${\cal C}$| introduced in subsection 3.5.

Because the radial computational domain is limited to r_min < r′ < r_max, the radial integration in equation (17) cannot be extended to r′ < r_min and r′ > r_max. However, neglecting the reprocessed starlight from outside the computational domain would result in an underestimate of the temperature near the computational boundaries. Near the inner computational boundary, an underestimated temperature amplifies the numerical wiggle in the radial disk structure (see appendix 3). To mitigate the underestimation of the temperature, we account for the flux outside the computational domain in an approximate way. Specifically, we write F_rep,↓ as

$$\begin{equation} F_{\rm rep, \downarrow } = F_{\rm rep,\! \downarrow ,\rm domain} + F_{\rm rep, \downarrow , in} + F_{\rm rep, \downarrow , out}, \end{equation}$$

(28)

where |$F_{\rm rep, \downarrow ,\rm domain}$| represents the flux from within the computational domain, whereas F_rep,↓,in and F_rep,↓,out represent the additional fluxes from r′ < r_min and r′ > r_max, respectively. Our prescription for these additional fluxes is described in appendix 3.

To suppress numerical instabilities arising from grid-scale fluctuations of the reprocessed starlight flux, we replace F_{rep, ↓}(r_j) at j = 2, 3, …, J − 1 with the geometric means of the raw fluxes at j − 1 and j + 1, |$\sqrt{F_{\rm rep, \downarrow }(r_{j-1})F_{\rm rep, \downarrow }(r_{j+1})}$|⁠. At j = 1 and J, we use |$\sqrt{F_{\rm rep, \downarrow }(r_{1})F_{\rm rep, \downarrow }(r_{2})}$| and |$\sqrt{F_{\rm rep, \downarrow }(r_{J-1})F_{\rm rep, \downarrow }(r_{J})}$| instead. We expect that this grid-scale smoothing would become unnecessary once we include the radial diffusion of the disk’s thermal emission in equation (21) in future work.

4 Validation with gapped disks

In this section, we test our global two-layer radiative transfer model against disks with an annular surface density gap. Specifically, we examine if the global two-layer model reproduces the temperature profiles of gapped disks previously obtained by Jang-Condell and Turner (2012) using a Monte Carlo radiative transfer model.

4.1 Model

Following Jang-Condell and Turner (2012), we consider axisymmetric gapped disks with the radial gas surface density profile given by

$$\begin{equation} \Sigma \!=\! 19 \biggl ({ \displaystyle \frac{r}{10\, \rm au}} \biggr ) ^{-0.9} \left\lbrace 1\!-\! d\exp \left[-\frac{(r\!-\! a)^2}{2w^2}\right]\right\rbrace \, \rm g\, cm^{-2}, \end{equation}$$

(29)

where |$19({r}/{10\, \rm au})^{-0.9}\, \rm g\, cm^{-2}$| is the surface density profile with no planet, |$a = 10\, \rm au$| is the planet’s orbital radius, and d and w represent the depth and width of the planet-carving gap, respectively. They considered one case with no planet and two cases with a planet of mass M_p = 70 or 200M_⊕ orbiting at 10 au from the central star. The parameters (d, w) were taken to be (0.56, 0.11a) and (0.84, 0.17a) for M_p = 70 and 200 M_⊕, respectively. The top row of figure 6 shows the gas surface density profiles for three disk models. The central star was assumed to have mass M_* = M_⊙, radius R_* = 2.6 R_⊙, surface temperature |$T_* = 4280\, \rm K$|⁠.

$Steady-state irradiation structure of the Jang-Condell and Turner (2012) disk models, reproduced with global two-layer radiative transfer calculations. Columns (a), (b), and (c) are for the disk models with no planet, a 70 M⊕ planet at 10 au, and a 200 M⊕ one. The planet is at 10 au. Top row: assumed gas surface density profiles (solid lines) compared with the no-planet profile (dashed). Second row: steady-state temperature profiles (thick solid lines) and equilibrium temperatures from Monte Carlo radiative transfer calculations (Jang-Condell & Turner 2012, red dashed lines). The thin solid lines give temperatures from the local two-layer model, equation (7), using μ* from the global two-layer model. Third row: fractional temperature difference between the global two-layer and Monte Carlo calculations. Thick and thin lines show the errors inside and outside the region r ≈ 7–17 au, where the Monte Carlo calculation is least affected by the computational domain boundaries. Fourth row: the heights of the starlight irradiation surface and infrared thermal photosphere, zs and zIR (thick and thin solid lines, respectively), normalized by r. The upper and lower green dotted lines indicate that zs ≈ 3–4 H. Bottom row: starlight grazing angle α* = sin −1μ*.$

Fig. 6.

Steady-state irradiation structure of the Jang-Condell and Turner (2012) disk models, reproduced with global two-layer radiative transfer calculations. Columns (a), (b), and (c) are for the disk models with no planet, a 70 M_⊕ planet at 10 au, and a 200 M_⊕ one. The planet is at 10 au. Top row: assumed gas surface density profiles (solid lines) compared with the no-planet profile (dashed). Second row: steady-state temperature profiles (thick solid lines) and equilibrium temperatures from Monte Carlo radiative transfer calculations (Jang-Condell & Turner 2012, red dashed lines). The thin solid lines give temperatures from the local two-layer model, equation (7), using μ_* from the global two-layer model. Third row: fractional temperature difference between the global two-layer and Monte Carlo calculations. Thick and thin lines show the errors inside and outside the region r ≈ 7–17 au, where the Monte Carlo calculation is least affected by the computational domain boundaries. Fourth row: the heights of the starlight irradiation surface and infrared thermal photosphere, z_s and z_IR (thick and thin solid lines, respectively), normalized by r. The upper and lower green dotted lines indicate that z_s ≈ 3–4 H. Bottom row: starlight grazing angle α_* = sin ⁻¹μ_*.

Open in new tab Download slide

Jang-Condell and Turner (2012) provided the temperature structure of the three disk models (M_p = 0, 70, and 200 M_⊕) using two radiative transfer models. One is a model conceptually similar to ours, relying on the locally one-dimensional analytic solution for a plane-parallel disk (Jang-Condell & Sasselov 2003, 2004; Jang-Condell 2008, 2009). The other is the Monte Carlo model developed by Turner et al. (2012), which directly follows the emission, absorption, and scattering of photon packets using the frequency-dependent opacities provided by Jang-Condell (2009). The two approaches returned similar results as shown in figure 4 of Jang-Condell and Turner (2012), so we only use their Monte Carlo results in the following comparison. The dashed lines in the second row of figure 6 show the temperature profiles from the Monte Carlo calculations by Jang-Condell and Turner (2012). We note that their Monte Carlo calculations only cover 3–20.8 au.

We adopt the mean opacities for stellar and disk thermal radiation computed from the frequency-dependent opacities of Jang-Condell (2009). For disk thermal radiation, we evaluate the mean opacities at a fixed temperature of |$T = 50\, \rm K$|⁠. The adopted mean opacities are |$\chi _* = 11.3\, \rm cm^{2}\, g^{-1}$|⁠, |$\kappa _* = 1.36\, \rm cm^{2}\, g^{-1}$|⁠, |$\chi _{\rm P} = 1.92\, \rm cm^{2}\, g^{-1}$|⁠, |$\chi _{\rm R} = 1.98\, \rm cm^{2}\, g^{-1}$|⁠, and , |$\kappa _{\rm P} = 0.975\, \rm cm^{2}\, g^{-1}$|⁠, which yield ε_* = 0.12, q = 5.7, f_↓ = 0.24 (see figure 2), and ε = 0.49. The disk model considered here is optically thick to its own thermal emission, with τ_R ∼ 10 at |$r \sim 10\, \rm au$|⁠. The thermal relaxation time at |$r\sim 10\, \rm au$| is |$t_{\rm th} \sim k_{\rm B}\Sigma /(m_{\rm g}\sigma _{\rm SB}T^3 \tau _{\rm R}) \sim 10^2\, \rm yr$|⁠.

Unlike the radiative equilibrium Monte Carlo calculations of Jang-Condell and Turner (2012), our global two-layer calculations are fully time-dependent. The timestep in our calculations is fixed to be 1 yr, which is shorter than the thermal relaxation time in the computational domain. The initial temperature profile is set to be |$T= 120(r/1\, \rm au)^{-3/7}\, \rm K$|⁠. This choice is arbitrary, but is not far from the steady-state temperature profile for the model with no planet.

Our two-layer calculations use a computational domain covering r_min = 1 au to r_max = 100 au, divided into 200 logarithmically spaced radial cells. Our computational domain is wider than in the Monte Carlo calculations by Jang-Condell and Turner (2012), and therefore our results are likely less affected by the computational boundaries and by any outer edge to the disk. We neglect radiation from the parent molecular cloud by setting T_ex = 0 in equation (21). The additional reprocessed starlight fluxes from outside the computational domain, F_{rep, ↓, in} and F_{rep, ↓, out}, are included.

4.2 Results

Figure 7 shows the time evolution of the temperature T for three Jang-Condell and Turner (2012) disk models obtained from our global two-layer calculations. In the M_p = 0 and 70 M_⊕ models, the temperature profile relaxes into a steady state. During relaxation, the temperature exhibits a damped oscillation with a period of ≈700 yr. This oscillation period is comparable to the thermal relaxation time of the disk, which is ∼200–400 yr for the model with no planet. This implies that the observed oscillation is not an artifact but reflects the system’s thermal relaxation. In the M_p = 200M_⊕ disk model, an oscillation of a similar period persists at a constant amplitude. The oscillation propagates toward the star, indicating that it is likely due to the thermal wave instability (see section 5 for more details about the long-term behavior of this instability). In the following, we describe the results for each disk model in more detail.

Fig. 7.

Evolution of the temperature structure in time-dependent radiative transfer calculations for three Jang-Condell and Turner (2012) disk models from figure 6. The temperature oscillates with a period similar to the disk’s thermal relaxation timescale, in a pattern that propatates toward the star, likely due to the thermal wave instability (subsection 4.2).

Open in new tab Download slide

4.2.1 M_p = 0 model

The steady-state structure of the M_p = 0 model at t = 3 × 10⁴ yr is shown in the left column of figure 6. The steady-state temperature profile follows a power law (see figure 8).⁴ This is reasonable because we assume a power-law surface density profile. At r ≈ 7–17 au, the steady-state temperature matches the radiative equilibrium temperature from the Monte Carlo calculation by Jang-Condell and Turner (2012) to within |$4\%$| (see the third row of figure 6). The excellent match can only be achieved when we account for the scattering of stellar radiation; using f_↓ = 0.5 [for ε_* = 1; see equation (2)] instead of the correct value of f_↓ = 0.24 would overestimate the temperature by |$\approx\! 40\%$|⁠. At r ≲ 7 au and r ≳ 17 au, the Monte Carlo result deviates from the power-law profile. This is likely because of the relatively narrow computational domain and disk outer edge adopted in the Monte Carlo calculation. In the following, comparisons between our results and those of Jang-Condell and Turner (2012) are only made at r ≈ 7–17 au.

Fig. 8.

Steady-state temperature profile from a global two-layer radiative transfer calculation for the Jang-Condell and Turner (2012) disk model with no planet (solid line) and a power-law fit T ∝ r^−0.48 (dotted line). The dashed line is the equilibrium temperature profile from the Monte Carlo calculation by Jang-Condell and Turner (2012).

Open in new tab Download slide

The irradiation surface lies at z_s ∼ 3–4H (the fourth row of figure 6), and the starlight grazing angle α_* is ≈0.02–0.03 over the entire disk (the bottom row of figure 6). The small drop in α_* at the inner computational boundary is an artifact and depends on how we treat |$\tau _{*,\rm in}$| and F_rep,↓,in. As discussed in appendix 3, our prescriptions for |$\tau _{*,\rm in}$| and F_rep,↓,in already suppress this inner boundary artifact significantly, if not completely. A similar inner boundary artifact can also be seen in the results for M_p = 70 and 200 M_⊕ [see columns (b) and (c) of figure 6].

4.2.2 M_p = 70M_⊕ model

The steady-state structure of the M_p = 70M_⊕ disk model at t = 3 × 10⁴ yr is shown in the center column of figure 6. The temperature profile has a dip and a bump around r ≈ 9 and 12 au, approximately corresponding to the bottom and outer edge of the planetary gap, respectively. As already pointed out by Jang-Condell and Turner (2012), these features arise because the gap’s trough is less exposed to stellar radiation and the gap’s outer rim is more illuminated. This can also be seen in the bottom row of figure 6, which shows that α_* has a local minimum and a maximum in the gap’s trough and outer rim, respectively.

At r ≈ 7–17 au, our temperature profile for the M_p = 70M_⊕ model matches the Monte Carlo result by Jang-Condell and Turner (2012) to an accuracy of |$\lesssim\!\! 15\%$|⁠. In particular, the error in the gap region falls below a few percents. The error is larger beyond the gap’s outer edge, reaching |$15\%$| at |$r \sim 15\, \rm au$|⁠.

4.2.3 M_p = 200M_⊕ model

Because the disk structure for the M_p = 200M_⊕ model is time-dependent due to the thermal wave instability, we generate steady-state disk structure by averaging the simulation result over t = 1.5–3 × 10⁴ yr. The obtained steady-state disk structure is shown in the right column of figure 6. In the inner half of the gap, the starlight grazing angle falls below 10⁻², indicating that this part almost falls into a shadow cast by the inner disk. As in the case of M_p = 70M_⊕, our two-layer model reproduces the gap temperature from the Monte Carlo calculation by Jang-Condell and Turner (2012) to an accuracy of a few %. The agreement is less good beyond the gap’s outer edge, where the error approaches |$20\%$|⁠.

It is interesting to ask why the M_p = 200 M_⊕ model is thermally unstable while the other two models are not. In general, the thermal wave instability operates in disks with large z_s/H (Wu & Lithwick 2021; Ueda et al. 2021). As shown in the upper panel of figure 9, the M_p = 200M_⊕ model is indeed the one among the three models that has the largest z_s/H, except around the gap’s sunlit outer edge (⁠|$r \sim 12\, \rm au$|⁠). In this model, the deep gap modifies the disk structure such that H/r becomes radially flat at r ∼ 8 au and r ∼ 20 au (see the lower panel of figure 9). As shown by Garaud and Lin (2007) and Wu and Lithwick (2021), a radially flat H/r generally leads to a high value of z_s/H. We suspect that the relatively high z_s/H around the 200M_⊕ planet’s gap destabilizes the disk.

Fig. 9.

Steady-state radial profiles of z_s/H and H/r (upper and lower panels) obtained using the global two-layer method for the three Jang-Condell and Turner (2012) disk models from figure 6. Note that high values of z_s/H are associated with plateaus of H/r (Garaud & Lin 2007; Wu & Lithwick 2021).

Open in new tab Download slide

4.2.4 The role of radial radiation transfer

We have shown that our global two-layer model reproduces the temperature profiles from the Monte Carlo calculations by Jang-Condell and Turner (2012) to within an accuracy of 20%. Before closing this section, we also show that local two-layer models would never produce such an accurate temperature profile for gapped disks. In the second row of figure 6, we overplot the temperature profiles that one would have from the local two-layer model, equation (7), using the profiles of μ_* derived from our global two-layer calculations. It can be seen that the local model appreciably underestimates and overestimates the temperatures at the gap’s trough and outer rim, respectively. This suggests that our global treatment of the reprocessed starlight is essential for accurately predicting the gap temperature.

5 Application to the thermal wave instability

The thermal wave instability is one of the most interesting targets of our global two-layer radiative transfer model. The previous simulations by Watanabe and Lin (2008), Ueda, Flock, and Birnstiel (2021), and Wu and Lithwick (2021) showed that the nonlinear stage of the instability is characterized by a train of inward moving temperature peaks with extended shadows. However, the simulations by Watanabe and Lin (2008) and Ueda, Flock, and Birnstiel (2021) used a simplified radiative transfer model with ad hoc radial smoothing for reprocessed starlight, which may not accurately resolve the thermal waves. Moreover, their simulations overestimated the disk cooling rate in optically thick regions as pointed out by Wu and Lithwick (2021). In contrast, Wu and Lithwick (2021) used a Monte Carlo radiative transfer code that fully includes the radial radiation transfer of reprocessed starlight. However, they only simulated the relatively early phase of the instability over several relaxation times.

In this section, we use our global two-layer model to simulate the long-term (t ≫ t_th) evolution of the thermal wave instability. The aims here are to study how the delayed thermal relaxation in optically thick regions, radial radiative transfer of reprocessed starlight, and radial numerical resolution affect the nonlinear development of the instability. We describe the model in subsection 5.1 and present the results in subsection 5.2.

5.1 Model

We adopt the disk model used by Ueda, Flock, and Birnstiel (2021). The gas surface density profile is given by |$\Sigma = 1700(r/1\, {\rm au})^{-3/2}\exp (-r/100\, \rm au)\, g\, cm^{-2}$|⁠. The disk opacities are assumed to scale with the dust-to-gas mass ratio f_d2g of opacity-dominating dust grains. For the fiducial value of f_d2g = 0.01, the mean opacities are |$\chi _{*} = \kappa _{*} = 8\, \rm cm^2\, g^{-1}$| and |$\chi _{\rm R} =\chi _{\rm P} = \kappa _{\rm P} = 4\, \rm cm^2\, g^{-1}$|⁠. For dust grains this is a free parameter of the model. In this study, we only consider this fiducial model. Because ε_* = 1, we have f_↓ = 0.5. External thermal radiation of |$T_{\rm ex} = 10\, \rm K$| is included, so that the disk temperature never falls below 10 K.

The radial computational domain ranges between 0.03 and 300 au and is divided into 480 logarithmically spaced cells, resulting in a radial resolution of dr/r = 0.019. The adopted radial resolution is two times better than in the simulation by Ueda, Flock, and Birnstiel (2021). In appendix 5, we also present simulation runs with lower resolution to study the convergence of our simulation results.

Our simulations differ from those by Ueda, Flock, and Birnstiel (2021) in the treatments of radial radiative transfer and thermal relaxation. To discriminate between the effects of the two different treatments, we also conduct a simulation using the same local radiative transfer model as in the simulations by Ueda, Flock, and Birnstiel (2021) but including the correction for thermal relaxation in optically thick regions, the factor (1 + 3τ_{R, mid}/4)⁻¹ in equation (22). Their radiative transfer model, originally developed by Watanabe and Lin (2008), approximates the downward reprocessed starlight flux as

$$\begin{equation} F_{\rm rep, \downarrow }(r) = \frac{L_*}{8\pi }\left\langle \frac{A_s(r^{\prime })}{r^{\prime 2}} + \frac{4R_*}{3\pi r^{\prime 3}}\right\rangle , \end{equation}$$

(30)

$$\begin{equation} A_s (r) = 1-\exp \left[-\int _{z_{\rm s}(r)}^\infty \kappa _{*}\rho (r,z^{\prime })dz^{\prime }\right], \end{equation}$$

(31)

where the angled brackets denote a radial average over a radial zone near r′ = r. This averaging is introduced to mimic the radial propagation of reprocessed starlight over distance ∼z_s(r). Following Ueda, Flock, and Birnstiel (2021), we use the Gaussian weight function exp [−(r′ − r)²/z_s(r)²] for the radial averaging in equation (30).

The simulations are carried out over 1 Myr. The simulated time covers ∼100t_th at 10 au and ∼10t_th at 1 au (see subsection 5.2). In comparison, the simulation reported by Wu and Lithwick (2021) only covers several thermal times at 10 au and less than one thermal time at 1 au. Therefore, our simulation for the first time fully captures thermal waves traveling from the 10 au to 1 au regions, with a correct treatment for thermal relaxation in the optically thick regime. The timestep Δt is taken to be 1 yr, which is 10 times shorter than the minimum thermal relaxation time in the disk (⁠|$\sim 10\, \rm yr$|⁠; see figure 11 in subsection 5.2). Using a shorter timestep of Δt = 0.25 yr makes no appreciable change in the simulation results.

Fig. 10.

Radial profiles of the disk interior temperature T (top panels), normalized irradiation surface height z_s/r (middle panels), and starlight grazing angle α_* (bottom panels) from the global two-layer calculation for the Ueda, Flock, and Birnstiel (2021) disk model with f_d2g = 0.01. The left and right columns show the time variation in the intervals t = 0.38–0.40 Myr and 0.4–1.0 Myr, respectively. The arrows in the top panels mark five selected temperature peaks (1, 2, 3, 1′, and 2′) migrating inward. At t = 0.40 Myr, peak 2 catches up with peak 1 and starts merging with it. The dotted lines in the middle panels indicate 4H/r.

Open in new tab Download slide

$Radial profiles of the thermal relaxation timescale tth [equation (27); solid line], local Keplerian time (dotted line), and vertical infrared optical depth to the midplane, τR, mid = κRΣ/2 (dashed line), for the Ueda, Flock, and Birnstiel (2021) disk model with fd2g = 0.01. The thermal relaxation time shown here uses the temperature distribution at $t = 1\, \rm Myr$ obtained from our simulation (see the top right panel of figure 10).$

Fig. 11.

Radial profiles of the thermal relaxation timescale t_th [equation (27); solid line], local Keplerian time (dotted line), and vertical infrared optical depth to the midplane, τ_{R, mid} = κ_RΣ/2 (dashed line), for the Ueda, Flock, and Birnstiel (2021) disk model with f_d2g = 0.01. The thermal relaxation time shown here uses the temperature distribution at |$t = 1\, \rm Myr$| obtained from our simulation (see the top right panel of figure 10).

Open in new tab Download slide

All our simulations assume vertical hydrostatic equilibrium (see subsection 3.4). As shown in subsection 5.2, the requirement t_th ≫ t_K for the hydrostatic equilibrium is fulfilled over the entire thermally unstable region.

5.2 Results

Figure 10 presents snapshots of the temperature distribution at selected times. We find that the temperature structure at r ≈ 0.3–30 au is unstable and exhibits oscillations (thermal waves) that propagate inward. The thermal waves consist of sharp temperature peaks with width ∼z_s(r_peak), where r_peak is the radial position of each peak. Each peak casts a shadow with a radial width ∼r_peak (see the middle and bottom panels of figure 10), and each shadow causes a dip in the temperature profile. The irradiation surface height in the unstable region exceeds 4H, which is higher than in the Jang-Condell and Turner (2012) disk models. This provides further support for the prediction by Wu and Lithwick (2021) that disks with larger z_s/H are more prone to the thermal wave instability.

At |$r \gtrsim 30\, \rm au$| and |$r \lesssim 0.3\, \rm au$|⁠, the thermal wave instability is suppressed for different reasons. In the outer region of |$r \gtrsim 30\, \rm au$|⁠, the external radiation flux |$\sigma _{\rm SB}T_{\rm ex}^4$| is comparable to or even dominates the reprocessed starlight flux F_{rep, ↓}, directly suppressing the thermal wave instability (Ueda et al. 2021). In the inner region of |$r \lesssim 0.1\, \rm au$|⁠, the finite size of the central star determines the starlight grazing angle [i.e., α_* ≈ 4R_*/(3πr); see equation (9)], which stabilizes the thermal wave instability as it is triggered by the variation of the stellar grazing angle with temperature fluctuations. In addition, the thermal timescale in this inner region is too long for thermal waves to develop within 1 Myr. Indeed, figure 11 shows that the thermal relaxation time t_th [equation (27)] at |$r \lesssim 0.1\, \rm au$| exceeds 1 Myr. Figure 11 also shows that the condition t_th ≫ t_K for vertical hydrostatic equilibrium is fulfilled in the entire region where the thermal wave instability is observed.

Compared to the thermal waves observed in the simulations by Ueda, Flock, and Birnstiel (2021; see their figures 2 and 3), our thermal waves are less coherent and propagate inward on different timescales at different radial positions. This can be more clearly seen in figure 12, which shows a spacetime plot of the temperature distribution. To highlight the wave components, we here normalize the temperature profile by the time-average T_avr(r) over t = 0.1–1 Myr. One can see that temperature peaks at smaller r migrate inward on longer timescales. The observed migration timescale at each r is crudely consistent with the local thermal timescale t_th shown in figure 11. This suggests that the radial variation of the migration speed is a manifestation of the radial variation of t_th. In contrast, in the simulations by Ueda, Flock, and Birnstiel (2021), t_th was nearly independent of r because they did not include the correction for the thermal relaxation rate. This explains why the thermal waves observed by Ueda, Flock, and Birnstiel (2021) propagate inward on a radially uniform timescale. For the Ueda, Flock, and Birnstiel (2021) disk model with f_d2g = 0.01, the thermal relaxation correction must be included because the disk is optically thick (⁠|$\tau _{\rm R,\rm mid} > 1$|⁠) at |$r \lesssim 100\, \rm au$| (see figure 11).

Fig. 12.

Spacetime plot for the temperature T, normalized by the time average T_avr(r) over t = 0.1–1 Myr, from the global two-layer calculation for the Ueda, Flock, and Birnstiel (2021) disk model with f_d2g = 0.01. This plot highlights the thermal waves.

Open in new tab Download slide

Another interesting feature that was not visible in the simulations by Ueda, Flock, and Birnstiel (2021) is the collision of temperature peaks. As an example, the top left panel of figure 10 shows how two temperature peaks labeled 1 and 2 collide. The outer peak 2 migrates faster than the inner peak 1, and they merge together at |$t \approx 0.40\, \rm Myr$|⁠. We speculate that the radially varying migration timescales of the temperature peaks induce their collisions.

Because the observed temperature peaks are narrow, it is important to examine whether the finite numerical resolution affects our simulation results. In appendix 5, we show that a radial resolution of dr/r ≲ 0.04 is enough to resolve the thermal waves.

Figure 13 shows a spacetime temperature plot from the simulation using the radiative transfer model of Watanabe and Lin (2008). The results are qualitatively similar to those from the global two-layer simulation presented above in that both feature collisions of temperature peaks. However, the Watanabe and Lin (2008) model produces wider temperature peaks than our global two-layer model. This is likely due to the ad hoc radial averaging of the reprocessed starlight flux adopted in the model. Therefore, we conclude that the radiative transfer model of Watanabe and Lin (2008) captures the qualitative nature of the thermal wave instability but should not be used for a quantitative study of the instability.

Fig. 13.

Same as figure 12, but from a simulation using the Watanabe and Lin (2008) radiative transfer model with Gaussian radial smoothing of the reprocessed starlight.

Open in new tab Download slide

5.3 Implications for dust evolution and planetesimal formation

As already noted by previous studies (Watanabe & Lin 2008; Ueda et al. 2021; Wu & Lithwick 2021), the thermal wave instability may have two important implications for dust evolution and planetesimal formation in dust-rich disks. Firstly, the thermal wave instability causes temporal variations in snow line locations. We demonstrate this in figure 14, where we show how the positions where T = 160, 70, and 20 K move with time in the simulation shown in figures 10 and 12. The selected temperatures correspond to the sublimation temperatures of H₂O, CO₂, and CO ices, respectively (Okuzumi et al. 2016). Overall, the plot indicates that the instability causes order-of-unity time variations in the positions of the snow lines. For instance, the H₂O snow line would migrate between ≈0.4–1.5 au on the timescale of ∼0.1 Myr. The previous simulation by Ueda, Flock, and Birnstiel (2021) predicted a migration timescale of ∼10 yr for the H₂O snow line. Our results indicate that the simulations by Ueda, Flock, and Birnstiel (2021) underestimated the migration timescale of the H₂O snow line because their simulations did not include the correction for thermal relaxation in the optically thick limit. Our new simulation suggests that the H₂O snow line can in fact move on a timescale comparable to the planet formation timescale, ∼0.1–1 Myr, if the disk’s optical thickness is large enough. We expect that the snow line migration should have important effects on planetesimal formation around the snow line (Ros & Johansen 2013; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Ida et al. 2021; Hyodo et al. 2021), in particular on the composition of forming planetesimals. In future work, we will explore the potential effects by including dust evolution in our global two-layer radiative transfer calculations.

Fig. 14.

Spacetime plot showing how the positions where T = 160 (thick lines), 70 (center thin lines), and 20 K (right thin lines) vary with time, obtained from the simulation for the Ueda, Flock, and Birnstiel (2021) disk model with f_d2g = 0.01. These positions correspond to the snow lines of H₂O, CO₂, and CO ices, respectively.

Open in new tab Download slide

Secondly, the thermal wave instability might produce pressure maxima that can trap dust particles. Dust particles in a gas disk are known to drift in the direction of increasing gas pressure (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977). Specifically, the radial drift velocity of a dust particle in a gas disk is given by

$$\begin{equation} v_{{\rm d},r} = \frac{\Omega _{\rm K} t_{\rm stop}}{1+(\Omega _{\rm K} t_{\rm stop})^2}\frac{c_{\rm s}^2}{v_{\rm K}}\frac{\partial \ln P}{\partial \ln r}, \end{equation}$$

(32)

where t_stop is the stopping time of the particle, v_K = rΩ_K is the local Keplerian velocity, and P = ρk_BT/m_g is the gas pressure. The radial drift velocity is proportional to the radial pressure gradient ∂ln P/∂ln r = ∂ln ρ/∂ln r + ∂ln T/∂ln r, with negative and positive radial gradients leading to inward and outward drift, respectively. A pressure maximum acts as a trap for radially drifting particles because their drift velocity converges there. Watanabe and Lin (2008) already suggested that the thermal wave instability can potentially produce local pressure maxima. Our new simulations relying on a more rigorous radiative transfer model confirm this possibility. Figure 15 shows some snapshots of the radial distribution of ∂ln P/∂ln r at the midplane from our global two-layer simulation presented in figure 10. We find that some of the temperature peaks observed in figure 10 (e.g., peak 2′ shown in the right panel of figure 10) indeed reverse the sign of the pressure gradient locally. Less pronounced temperature peaks that do not reverse the pressure gradient may also slow down particle inward drift and thereby promote planetesimal formation. However, a question remains as to whether the steep pressure variations observed here are hydrodynamically stable, because they may violate the Rayleigh criterion for stable disk rotation (e.g., Yang & Menou 2010). More fundamentally, the actual thermal waves around the midplane of optically thick disks could be less pronounced than predicted from our calculations owing to the vertically varying thermal relaxation time (Pavlyuchenkov et al. 2022). Addressing these issues is beyond the scope of this paper but should be done in future work.

Fig. 15.

Radial pressure gradient ∂ln P/∂ ln r at the midplane as a function of r at different times, obtained from the global two-layer simulation presented in figure 10. The arrows indicate the positions of temperature peak 2′ shown in the top right panel of figure 10.

Open in new tab Download slide

6 Summary

We have presented a new two-layer radiative transfer model for computing the radial temperature distribution of axisymmetric, irradiated protoplanetary disks. Unlike the standard two-layer model, our new model explicitly accounts for the radial transfer of reprocessed starlight and is therefore applicable to disks with shadowed regions. Our global two-layer model is conceptually similar to the radiative transfer models of Jang-Condell and Sasselov (2003, 2004) and Jang-Condell (2008), but is much simpler and computationally more efficient thanks to the assumed disk axisymmetry. In addition, our model treats the radial temperature distribution as time-dependent (subsection 3.5), allowing us to study disks’ thermal instabilities.

We have tested the global two-layer model against the Monte Carlo radiative transfer calculations by Jang-Condell and Turner (2012) for gapped disks (section 4). We have found that the steady-state temperature profiles from our global two-layer calculations match the previous Monte Carlo results to within an accuracy of 20% even with a deep gap carved by a 200M_⊕ planet (figure 6). We have also found that disks with larger z_s/H are more prone to the thermal wave instability, confirming the prediction from the linear analysis by Wu and Lithwick (2021).

We have then used the global two-layer model to simulate, for the first time, the long-term evolution of the thermal wave instability with a correct treatment for thermal relaxation in the optically thick limit (section 5). Our simulations reveal that the thermal waves are highly chaotic, with inward migrating temperature peaks colliding and merging frequently (figures 10 and 12). Our simulations also show that the migration timescale of the temperature peaks increases as they move inward, reaching ≳0.1 Myr at 1 au. These observed properties of the thermal waves are likely attributed to the radially varying thermal relaxation timescale in our simulations. Snow lines also migrate on a timescale similar to the local thermal relaxation time (figure 14). Sharp temperature peaks produced by the thermal wave instability can reverse the sign of the radial pressure gradient locally (figure 15), indicating that they may act as a dust trap. Our future modeling will include dust evolution and study the coupled evolution of dust and disk temperature structure in detail.

Finally, we note that our current treatment of the disk vertical density structure is greatly simplified, relying on the isothermal and hydrostatic approximations (subsection 3.4). The isothermal approximation may introduce some errors in the evaluation of the irradiation surface height, although the errors appear to be small (Watanabe & Lin 2008). The hydrostatic approximation is inapplicable to the outer disk region where the thermal relaxation time is shorter than the orbital time. Adopting a two-temperature vertical density profile (Watanabe & Lin 2008) and treating the disk scale height as a dynamical variable (Dullemond 2000) may allow us to relax these assumptions without having to fully solve the hydrodynamic equation of motion. We will pursue this direction in future work. Including accretion heating, radial diffusion of the disk’s own thermal radiation, and vertically varying thermal relaxation timescales will be other important future directions.

Acknowledgements

The authors thank Tilman Birnstiel, Cornelis Dullemond, and Mario Flock for helpful discussions and Ralf Siebenmorgen for a careful review of the manuscript. This work was supported by JSPS KAKENHI Grant Numbers JP18H05438, JP19J01929, JP19K03926, JP20H00182, and JP20H01948. This work was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA and with the support of NASA grant 18-2XRP18_2-0059.

Appendix 1. Fraction of starlight reprocessed downward

Here, we derive the fraction of the stellar radiation reprocessed downward, f_↓ [equation (2)], using the locally plane-parallel disk model of Calvet et al. (1991).

We consider stellar (visible) radiation incident on a plane-parallel atmosphere at angle α_* = sin ⁻¹μ_* with respect to the atmosphere’s surface. Particles in the atmosphere either absorb or scatter the incident light and reprocess the absorbed component into infrared radiation. Part of the multiple-scattered starlight escapes from the atmosphere and the rest is eventually converted into infrared radiation. Calvet et al. (1991) derived the mean intensities and vertical fluxes of the multiple-scattered starlight and disk thermal emission in a locally plane-parallel disk in radiative equilibrium. They used the first and second moments of the frequency-integrated radiative transfer equations for the visual and infrared radiation with the Eddington approximation. Jang-Condell and Sasselov (2004) extended the model of Calvet et al. (1991) so that the model can deal with infrared scattering. The moment equations for the scattered starlight read [see also equations (10) and (11) of D’Alessio et al. 1998

$$\begin{equation} \frac{dJ_{\rm s}}{d\tau _{\rm s}} = \frac{3}{4\pi } F_{\rm s}, \end{equation}$$

(A1)

$$\begin{equation} \frac{dF_{\rm s}}{d\tau _{\rm s}} = 4\pi \epsilon _* J_{\rm s} - (1-\epsilon _*) F_*e^{-\tau _{\rm s}/\mu _*}, \end{equation}$$

(A2)

where J_s and F_s are the frequency-integrated mean intensity and vertical upward flux for the scattered starlight, respectively, |$\tau _{\rm s} = \int _z^\infty \chi _* \rho dz^{\prime }$| is the vertical extinction optical depth, and ε_* is the ratio between the absorption and extinction opacities κ_* and χ_* for the stellar radiation [see equation (3)]. We follow D’Alessio et al. (1998) and evaluate κ_* and χ_* using the Planck means at the stellar surface temperature T_*,

$$\begin{equation} \kappa _* = \frac{\pi }{\sigma _{\rm SB} T_*^4} \int _0^\infty \kappa _\nu B_\nu (T_*) d\nu , \end{equation}$$

(A3)

$$\begin{equation} \chi _* = \frac{\pi }{\sigma _{\rm SB} T_*^4} \int _0^\infty \chi _\nu B_\nu (T_*) d\nu , \end{equation}$$

(A4)

where κ_ν and χ_ν are the monochromatic absorption and extinction opacities at frequency ν, respectively, and B_ν is the Planck function.

The moment equations for the disk’s infrared radiation are [see also equations (15) and (16) of D’Alessio et al. 1998]

$$\begin{equation} \frac{dJ_{\rm d}}{d\tau } = \frac{3}{4\pi }F_{\rm d}, \end{equation}$$

(A5)

$$\begin{equation} \frac{dF_{\rm d}}{d\tau } = 4\pi \epsilon (J_{\rm d} - B_{\rm d}), \end{equation}$$

(A6)

where J_d, F_d, and B_d are the frequency-integrated mean intensity, vertical upward flux, and Planck function for the disk thermal radiation, respectively, |$\tau = \int _z^\infty \chi \rho dz^{\prime }$| is the mean vertical optical depth for the infrared radiation, and ε = κ/χ. The mean opacities κ and χ appearing here are formally defined as the intensity- and flux-weighted averages of the absorption and extinction opacities for the disk thermal radiation, respectively. Following D’Alessio et al. (1998), we approximate κ and χ with the Planck and Rosseland mean opacities,

$$\begin{equation} \kappa _{\rm P} = \frac{\pi }{\sigma _{\rm SB} T^4} \int _0^\infty \kappa _\nu B_\nu (T) d\nu , \end{equation}$$

(A7)

$$\begin{equation} \chi _{\rm R} = \frac{4\sigma _{\rm SB} T^3}{\pi \int _0^\infty \chi _\nu ^{-1} {[dB_\nu (T)}/{dT}] d\nu } \end{equation}$$

(A8)

(see Hubeny et al. 2003 for the validity of these choices at large optical depths).

To fix the boundary conditions, we assume that both F_s and F_d vanish at large optical depths and that the outgoing radiation at the disk surface is hemispherically isotropic:

$$\begin{equation} F_{\rm s}(\tau _{\rm s}=0) = 2\pi J_{\rm s}(\tau _{\rm s}=0), \end{equation}$$

(A9)

$$\begin{equation} F_{\rm d}(\tau =0) = 2\pi J_{\rm d}(\tau =0). \end{equation}$$

(A10)

The resulting J_d(τ_s) is given by equation (12) of Calvet et al. (1991). For μ_* ≪ 1 and |$\tau _{\rm eff} \equiv \sqrt{3\epsilon } \tau \gg 1$|⁠, the expression for J_d can be greatly simplified as

$$\begin{equation} J_{\rm d} \approx \frac{1}{2\pi } \left[\epsilon _* + \left(2\sqrt{3\epsilon _*}+\frac{3}{q}\right)\frac{1-\epsilon _*}{3+2\sqrt{3\epsilon _*}} \right] \mu _* F_*. \end{equation}$$

(A11)

Because the thermal radiation at these depths is nearly isotropic, we may decompose it into upward and downward components that are hemispherically isotropic. The downward component has a flux of magnitude F_{rep, ↓} = πJ_d, which can be rewritten as F_{rep, ↓} = f_↓μ_*F_* if we define f_↓ by equation (2).

Appendix 2. Reprocessed starlight flux for plane-parallel disks

In this section, we prove that our expression for F_rep,↓ [equation (18)] recovers the more familiar expression F_rep,↓ = f_↓μ_*F_* [equation (1)] in the limit where h is small compared to the other relevant radial length scales. In this limit, the Lorentzian factor 1/[(r′ − r)² + h²] in the weighting function W is sharply peaked at r′ = r, and therefore the other factors inside the radial integration in equation (18) can be approximately evaluated at r′ = r. We thus obtain

$$\begin{eqnarray} F_{\rm rep, \downarrow }(r) \approx \frac{ f_\downarrow (r) \mu _*(r)F_*(r) h_{r=r^{\prime }} }{\pi } \int _0^\infty \frac{dr^{\prime }}{(r^{\prime }-r)^2 +h_{r=r^{\prime }}^2} ,\nonumber\\ \end{eqnarray}$$

(A12)

where we have used that h ≪ r ≈ r′, E(1) ≈ 1, and |$\sqrt{1+(dz_{\rm s}/dr^{\prime })^2} \approx 1$|⁠. The radial integration in equation (A12) yields

$$\begin{eqnarray} \int _0^\infty \frac{dr^{\prime }}{(r^{\prime }-r)^2 +h_{r=r^{\prime }}^2} = \frac{1}{2h_{r=r^{\prime }} }\left[\pi + 2\tan ^{-1} \biggl ({ \displaystyle \frac{r}{h_{r^{\prime } =r}}} \biggr ) \right].\nonumber\\ \end{eqnarray}$$

(A13)

For |$h_{r^{\prime } =r} \ll r$|⁠, we have |$\tan ^{-1}(r/h_{r^{\prime } =r}) \approx \pi /2$|⁠, which yields F_{rep, ↓} ≈ f_↓μ_*F_* .

Appendix 3. Numerical treatment of the computational boundaries

As noted in section 3, our radiative transfer calculations require the starlight optical depth to the inner computational boundary, |$\tau _{*,\rm in}$|⁠, and the reprocessed starlight fluxes from outside the inner and outer computational boundaries, F_{rep, ↓, in} and F_{rep, ↓, out}. All these depend on the conditions outside the computational domain and therefore must be assumed.

For |$\tau _{*,\rm in}$|⁠, we adopt the following prescription:

$$\begin{equation} \tau _{*,\rm in} = A_\tau \chi _* \rho (r_{\rm min},r_{\rm min} \tan \theta ) \frac{r_{\rm min}}{\cos \theta }. \end{equation}$$

(A14)

Here, r_min/cos θ stands for the distance from the coordinate origin to the position (r, rtan θ). The dimensionless number A_τ is a free parameter and is expected to be of order unity for uniformly flared disks with smooth radial density structure. Flock et al. (2016) adopted a similar prescription for τ_{*, in} with A_τ = 1. However, as shown below, the choice A_τ = 1 introduces a small wiggle in the temperature distribution near the inner computational boundary for the simulations presented in section 4. Throughout this paper, we adopt A_τ = 0.3, which better suppresses the wiggle for these particular simulations.

For F_{rep, ↓, in} and F_{rep, ↓, out}, we approximate the irradiation surface at r′ < r_min and r > r_max with flat planes lying at heights z_s(r_min) and z_s(r_max), extending radially over 0 < r′ < r_min and r′ > r_max, and producing vertical fluxes of the reprocessed starlight per irradiation surface area of |$[f_\downarrow \mu _* F_*]_{r^{\prime }=r_{\rm min}}$| and |$[f_\downarrow \mu _* F_*]_{r^{\prime }=r_{\rm max}}$| respectively. With these simplifications, the additional fluxes can be analytically calculated as [see equations (20) and (21) of Ueda, Okuzumi, and Flock 2017 for a derivation]

$$\begin{eqnarray} F_{\rm rep, \downarrow , in}(r) &=& \frac{\left[f_\downarrow \mu _* F_*\right]_{r^{\prime }=r_{\rm min}}}{\pi } \int _0^{r_{\rm min}} \frac{h_{r^{\prime }=r_{\rm min}}dr^{\prime }}{h_{r^{\prime }=r_{\rm min}}^2+(r-r^{\prime })^2} \nonumber \\ &=& \frac{\left[f_\downarrow \mu _* F_*\right]_{r^{\prime }=r_{\rm min}}}{\pi } \bigg[\tan ^{-1} \biggl ({ \displaystyle \frac{r}{h_{r^{\prime } = r_{\rm min}}}} \biggr ) \nonumber\\ &&-\tan ^{-1} \biggl ({ \displaystyle \frac{r-r_{\rm min}}{h_{r^{\prime } = r_{\rm min}}}} \biggr ) \bigg], \end{eqnarray}$$

(A15)

$$\begin{eqnarray} F_{\rm rep, \downarrow , out}(r) &=& \frac{[f_\downarrow \mu _* F_*]_{r^{\prime }=r_{\rm max}}}{\pi } \int _{r_{\rm max}}^\infty \frac{h_{r^{\prime }=r_{\rm max}}dr^{\prime }}{h_{r^{\prime }=r_{\rm max}}^2+(r-r^{\prime })^2} \nonumber \\ &=& \frac{[f_\downarrow \mu _* F_*]_{r^{\prime }=r_{\rm max}}}{\pi } \left[\frac{\pi }{2}-\tan ^{-1} \biggl ({ \displaystyle \frac{r_{\rm max}-r}{h_{r^{\prime } = r_{\rm max}}}} \biggr ) \right].\nonumber\\ \end{eqnarray}$$

(A16)

Below, we examine how our prescriptions for |$\tau _{*,\rm in}$|⁠, F_rep,↓,in, and F_rep,↓,out affect the calculated disk structure near the boundaries. We select the Jang-Condell and Turner (2012) disk model with no planet for a case study. The left column of figure 16 shows the steady-state structure near the inner boundary of this disk model calculated with different values of |$A_{\tau ,\rm in}$|⁠. One can see that the choices |$A_{\tau , \rm in} = 0.1$| and 1 introduce a wiggle in the radial temperature and grazing angle profiles. This artifact arises when a too large or small |$\tau _{*, \rm in}$| causes an underestimated or overestimated irradiation surface height z_s, respectively, near the inner boundaries. With our default choice |$A_{\tau , \rm in} = 0.3$| the wiggle is greatly suppressed, and the temperature profile follows a single power law T ∝ r^−0.48 down to the very vicinity of the inner boundary (see also figure 8).

$Steady-state disk structure of the Jang-Condell and Turner (2012) disk model with no planet obtained from global two-layer calculations with different prescriptions for the inner boundary. The upper and lower panels zoom in on the interior temperature T and starlight grazing angle α*, respectively, near the inner computational boundary. The left column compares the calculations with three different values of the parameter $A_{\tau , \rm in}$ controlling the starlight optical depth to the inner boundary [see equation (A14)], with $A_{\tau , \rm in} = 0.3$ being the default value. The right column compares the calculations for $A_{\tau , \rm in} = 0.3$ with and without the additional reprocessed flux beyond the inner boundary, Frep,↓,in [equation (A15)]. The dotted line in the upper panels is a power-law fit for the bulk temperature profile in the computational domain shown in figure 8.$

Fig. 16.

Steady-state disk structure of the Jang-Condell and Turner (2012) disk model with no planet obtained from global two-layer calculations with different prescriptions for the inner boundary. The upper and lower panels zoom in on the interior temperature T and starlight grazing angle α_*, respectively, near the inner computational boundary. The left column compares the calculations with three different values of the parameter |$A_{\tau , \rm in}$| controlling the starlight optical depth to the inner boundary [see equation (A14)], with |$A_{\tau , \rm in} = 0.3$| being the default value. The right column compares the calculations for |$A_{\tau , \rm in} = 0.3$| with and without the additional reprocessed flux beyond the inner boundary, F_rep,↓,in [equation (A15)]. The dotted line in the upper panels is a power-law fit for the bulk temperature profile in the computational domain shown in figure 8.

Open in new tab Download slide

The additional reprocessed flux F_rep,↓,in significantly contributes to removing the inner boundary artifact. Without F_rep,↓,in, the wiggle near the inner boundary survives even with |$A_{\tau , \rm in} = 0.3$| as shown in the right column of figure 16.

Figure 17 shows how F_rep,↓,out works near the outer computational boundary. Without this additional flux, the temperature profile deviates downward from the power law T ∝ r^−0.48 by up to 30% toward the outer boundary. Our prescription for F_rep,↓,out reduces this deviation to |$\lesssim 1\%$|⁠. Therefore, unless the outer computational boundary represents the disk’s true outer edge, one should apply this prescription to prevent temperature underestimation.

Fig. 17.

Steady-state temperature distribution in the outer region of the Jang-Condell and Turner (2012) disk with no planet from global two-layer calculations with and without the additional reprocessed flux beyond the out boundary, F_rep,↓,out [equation (A16)]. The dotted line in the temperature plot is a power-law fit for the bulk temperature profile in the computational domain shown in figure 8.

Open in new tab Download slide

Appendix 4. Emissivity of an isothermal slab

Here, we present an approximate model for the emissivity and absorptivity of a disk using a slab of uniform temperature T. From Kirchhoff’s law for thermal radiation, the emissivity and absorptivity of this isothermal slab must be equal. This allows us to derive its emissivity by considering a special case where there is no incoming radiation to the slab.

To obtain the frequency-integrated mean intensity J_d and flux F_d for the slab’s thermal radiation, we use the moment equations (A5) and (A6) already presented in appendix 1, but here relax the assumption that the disk’s optical thickness is infinitely large. Assuming that there is no incoming radiation and the outgoing radiation is isotropic at the slab’s upper and lower boundaries, the radiation field must obey the boundary conditions

$$\begin{equation} F_{\rm d}(\tau =0) = 2\pi J_{\rm d}(\tau =0), \end{equation}$$

(A17)

$$\begin{equation} F_{\rm d}(\tau =2\tau _{\rm mid}) = -2\pi J_{\rm d}(\tau =2\tau _{\rm mid}) , \end{equation}$$

(A18)

where τ_mid is the extinction optical depth to the midplane and 2τ_mid stands for the extinction optical thickness of the whole slab. By solving the moment equations with the boundary conditions, we obtain

$$\begin{eqnarray} J_{\rm d} (\tau ) = B_{\rm d}\left[1 -\frac{\sqrt{3}\cosh (\tau _{\rm eff}-\tau _{\rm eff,mid})}{\sqrt{3}\cosh (\tau _{\rm eff,mid}) +2\sqrt{\epsilon } \sinh (\tau _{\rm eff,mid})}\right]\nonumber\\ \end{eqnarray}$$

(A19)

with |$\tau _{\rm eff} = \sqrt{3\epsilon }\tau$| and |$\tau _{\rm eff,mid} = \sqrt{3\epsilon }\tau _{\rm mid}$|⁠. A similar expression for the mean intensity was also derived by Miyake and Nakagawa (1993), Jang-Condell (2008), Inoue, Oka, and Nakamoto (2009), and Birnstiel et al. (2018) but their expression is slightly different from ours because they used the two-stream approximation for the boundary condition. We have used the hemispherically isotropic outgoing boundary condition [equations (A17) and (A18)] to be consistent with the plane-parallel disk models of Calvet et al. (1991) and Jang-Condell and Sasselov (2004) that we used in our derivation of f_↓ (appendix 1).

The thermal flux at the upper boundary is

$$\begin{equation} F_{\rm d} (\tau =0) = \frac{4\pi B_{\rm d}\sqrt{\epsilon }\tanh (\tau _{\rm eff,mid})}{\sqrt{3} +2\sqrt{\epsilon } \tanh (\tau _{\rm eff,mid})}. \end{equation}$$

(A20)

The emissivity, the ratio between F_d(τ = 0) and blackbody radiation flux πB_d, is

$$\begin{equation} \frac{F_{\rm d} (\tau =0)}{\pi B_{\rm d}} = \frac{4 \sqrt{\epsilon }\tanh (\tau _{\rm eff,mid})}{\sqrt{3} +2\sqrt{\epsilon } \tanh (\tau _{\rm eff,mid})}. \end{equation}$$

(A21)

Equation (A21) has limiting expressions:

$$\begin{eqnarray} { \displaystyle \frac{F_{\rm d} (\tau =0)}{\pi B_{\rm d}}} &\approx & \left\lbrace \begin{array}{ll}4\sqrt{{ \displaystyle \frac{\epsilon }{3}} } \tau _{\rm eff,mid}, & \quad \tau _{\rm eff,mid} \ll 1 ,\\ { \displaystyle \frac{4 \sqrt{\epsilon }}{\sqrt{3} +2\sqrt{\epsilon } }} , \quad & \quad \tau _{\rm eff,mid} \gg 1 . \end{array} \right. \end{eqnarray}$$

(A22)

Appendix 5. The thermal wave instability: radial resolution dependence

To examine whether the finite numerical resolution affects our simulation results for the thermal wave instability (section 5), we performed two additional simulation runs for the same disk model but with two and four times coarser resolutions of dr/r = 0.038 and dr/r = 0.077.

Figure 18 compares the spacetime plots of T/T₀ from the low-resolution runs with the corresponding plot from the original simulation with dr/r = 0.019 already shown in subsection 5.2 (figure 12). For dr/r = 0.077, small-scale oscillations are significantly suppressed and only become prominent after |$t \gtrsim 0.8\, \rm Myr$| at |$r \gtrsim 2\, \rm au$|⁠. No collisions of temperature peaks are observed, and the wave pattern is much more coherent than in higher-resolution runs. The temperature peaks are wider than in the original run, suggesting that even the large-scale oscillations are considerably affected by the low resolution in this run. In the intermediate-resolution run with dr/r = 0.038, small-scale fluctuations become visible from |$t \sim 0.2\, \rm Myr$| and the wave pattern after this time is similar to that seen in the original run. From this convergence study, we can conclude that our simulation is well converged at a resolution of dr/r ≲ 0.04.

Fig. 18.

Same as figure 12, but comparing the fiducial simulation with the radial resolution of dr/r = 0.019 (left panel, already shown in figure 12) with those with two and four times coarser resolutions dr/r = 0.038 and 0.077 (center and right panels).

Open in new tab Download slide

Footnotes

1

However, it is under debate whether silicates or water ice are stickier (e.g., Kimura et al. 2015; Gundlach & Blum 2015; Musiolik & Wurm 2019).

2

We avoid calling this the irradiation instability (Wu & Lithwick 2021) because this term was previously used for a different disk instability (Fung & Artymowicz 2014).

3

The quantity h can become negative because it compares z_s and z_IR at different radial locations. For instance, z_s at small r′ can be smaller than z_IR at large r.

4

Our best-fitting T ∝ r^−0.48 (see figure 8) is steeper than the well-known profile T ∝ r^−3/7 for optically thick disks with radially constant z_s/H (Chiang & Goldreich 1997). In our temperature calculations, the ratio z_s/H generally varies with r because z_s is computed from the radial optical depth. In this no-planet disk model, z_s decreases from ∼4H at |$r \approx 1\, \rm au$| to ∼3H at |$r \approx 25\, \rm au$| (see the bottom left panel of figure 6) and therefore the irradiation surface flares more slowly with r than in constant z_s/H models, explaining the steeper temperature profile.

References

Adachi

I.

,

Hayashi

C.

,

Nakazawa

K.

1976

,

Prog. Theor. Phys.

,

56

,

1756

Crossref

Search ADS

Akiyama

E.

,

Momose

M.

,

Hayashi

H.

,

Kitamura

Y.

2011

,

PASJ

,

63

,

1059

Crossref

Search ADS

Armitage

P. J.

2010

,

Astrophysics of Planet Formation

(

Cambridge

:

Cambridge University Press

)

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Banzatti

A.

,

Pinilla

P.

,

Ricci

L.

,

Pontoppidan

K. M.

,

Birnstiel

T.

,

Ciesla

F.

2015

,

ApJ

,

815

,

L15

Crossref

Search ADS

Béthune

W.

,

Latter

H.

2020

,

MNRAS

,

494

,

6103

Crossref

Search ADS

Birnstiel

T.

et al.

2018

,

ApJ

,

869

,

L45

Crossref

Search ADS

Bitsch

B.

,

Crida

A.

,

Morbidelli

A.

,

Kley

W.

,

Dobbs-Dixon

I.

2013

,

A&A

,

549

,

A124

Crossref

Search ADS

Bitsch

B.

,

Johansen

A.

,

Lambrechts

M.

,

Morbidelli

A.

2015

,

A&A

,

575

,

A28

Crossref

Search ADS

Calahan

J. K.

et al.

2021

,

ApJS

,

257

,

17

Crossref

Search ADS

Calvet

N.

,

Patino

A.

,

Magris

G. C.

,

D’Alessio

P.

1991

,

ApJ

,

380

,

617

Crossref

Search ADS

Cannizzo

J. K.

1993

,

ApJ

,

419

,

318

Crossref

Search ADS

Carrasco-González

C.

et al.

2019

,

ApJ

,

883

,

71

Crossref

Search ADS

Chiang

E. I.

,

Goldreich

P.

1997

,

ApJ

,

490

,

368

Crossref

Search ADS

Chiang

E. I.

,

Joung

M. K.

,

Creech-Eakman

M. J.

,

Qi

C.

,

Kessler

J. E.

,

Blake

G. A.

,

van Dishoeck

E. F.

2001

,

ApJ

,

547

,

1077

Crossref

Search ADS

Chokshi

A.

,

Tielens

A. G. G. M.

,

Hollenbach

D.

1993

,

ApJ

,

407

,

806

Crossref

Search ADS

Cuzzi

J. N.

,

Zahnle

K. J.

2004

,

ApJ

,

614

,

490

Crossref

Search ADS

D’Alessio

P.

,

Cantö

J.

,

Calvet

N.

,

Lizano

S.

1998

,

ApJ

,

500

,

411

Crossref

Search ADS

D’Alessio

P.

,

Calvet

N.

,

Hartmann

L.

2001

,

ApJ

,

553

,

321

Crossref

Search ADS

D’Alessio

P.

,

Calvet

N.

,

Hartmann

L.

,

Lizano

S.

,

Cantó

J.

1999

,

ApJ

,

527

,

893

Crossref

Search ADS

Davis

S. S.

2005

,

ApJ

,

620

,

994

Crossref

Search ADS

Dominik

C.

,

Tielens

A. G. G. M.

1997

,

ApJ

,

480

,

647

Crossref

Search ADS

Drążkowska

J.

,

Alibert

Y.

2017

,

A&A

,

608

,

A92

Crossref

Search ADS

Dullemond

C. P.

2000

,

A&A

,

361

,

L17

Dullemond

C. P.

,

Dominik

C.

2004a

,

A&A

,

417

,

159

Crossref

Search ADS

Dullemond

C. P.

,

Dominik

C.

2004b

,

A&A

,

421

,

1075

Crossref

Search ADS

Dullemond

C. P.

,

Dominik

C.

,

Natta

A.

2001

,

ApJ

,

560

,

957

Crossref

Search ADS

Dullemond

C. P.

,

Isella

A.

,

Andrews

S. M.

,

Skobleva

I.

,

Dzyurkevich

N.

2020

,

A&A

,

633

,

A137

Crossref

Search ADS

Dullemond

C. P.

,

Juhasz

A.

,

Pohl

A.

,

Sereshti

F.

,

Shetty

R.

,

Peters

T.

,

Commercon

B.

,

Flock

M.

2012

,

Astrophysics Source Code Library, ascl:1202.015

Dullemond

C. P.

,

van Zadelhoff

G. J.

,

Natta

A.

2002

,

A&A

,

389

,

464

Crossref

Search ADS

Flaherty

K. M.

,

Hughes

A. M.

,

Teague

R.

,

Simon

J. B.

,

Andrews

S. M.

,

Wilner

D. J.

2018

,

ApJ

,

856

,

117

Crossref

Search ADS

Flaherty

K.

et al.

2020

,

ApJ

,

895

,

109

Crossref

Search ADS

Flock

M.

,

Fromang

S.

,

Turner

N. J.

,

Benisty

M.

2016

,

ApJ

,

827

,

144

Crossref

Search ADS

Fung

J.

,

Artymowicz

P.

2014

,

ApJ

,

790

,

78

Crossref

Search ADS

Garaud

P.

,

Lin

D. N. C.

2007

,

ApJ

,

654

,

606

Crossref

Search ADS

Gundlach

B.

,

Blum

J.

2015

,

ApJ

,

798

,

34

Crossref

Search ADS

Hirose

S.

,

Turner

N. J.

2011

,

ApJ

,

732

,

L30

Crossref

Search ADS

Hubeny

I.

,

Burrows

A.

,

Sudarsky

D.

2003

,

ApJ

,

594

,

1011

Crossref

Search ADS

Hyodo

R.

,

Guillot

T.

,

Ida

S.

,

Okuzumi

S.

,

Youdin

A. N.

2021

,

A&A

,

646

,

A14

Crossref

Search ADS

Hyodo

R.

,

Ida

S.

,

Charnoz

S.

2019

,

A&A

,

629

,

A90

Crossref

Search ADS

Ida

S.

,

Guillot

T.

2016

,

A&A

,

596

,

L3

Crossref

Search ADS

Ida

S.

,

Guillot

T.

,

Hyodo

R.

,

Okuzumi

S.

,

Youdin

A. N.

2021

,

A&A

,

646

,

A13

Crossref

Search ADS

Inoue

A. K.

,

Oka

A.

,

Nakamoto

T.

2009

,

MNRAS

,

393

,

1377

Crossref

Search ADS

Isella

A.

,

Turner

N. J.

2018

,

ApJ

,

860

,

27

Crossref

Search ADS

Jang-Condell

H.

2008

,

ApJ

,

679

,

797

Crossref

Search ADS

Jang-Condell

H.

2009

,

ApJ

,

700

,

820

Crossref

Search ADS

Jang-Condell

H.

,

Sasselov

D. D.

2003

,

ApJ

,

593

,

1116

Crossref

Search ADS

Jang-Condell

H.

,

Sasselov

D. D.

2004

,

ApJ

,

608

,

497

Crossref

Search ADS

Jang-Condell

H.

,

Turner

N. J.

2012

,

ApJ

,

749

,

153

Crossref

Search ADS

Kenyon

S. J.

,

Hartmann

L.

1987

,

ApJ

,

323

,

714

Crossref

Search ADS

Kim

S.

,

Nomura

H.

,

Tsukagoshi

T.

,

Kawabe

R.

,

Muto

T.

2019

,

ApJ

,

872

,

179

Crossref

Search ADS

Kimura

H.

,

Wada

K.

,

Senshu

H.

,

Kobayashi

H.

2015

,

ApJ

,

812

,

67

Crossref

Search ADS

Krügel

E.

2008

,

An Introduction to the Physics of Interstellar Dust

(

New York

:

Taylor & Francis

)

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Kusaka

T.

,

Nakano

T.

,

Hayashi

C.

1970

,

Prog. Theor. Phys.

,

44

,

1580

Crossref

Search ADS

Latter

H. N.

,

Balbus

S.

2012

,

MNRAS

,

424

,

1977

Crossref

Search ADS

Lynden-Bell

D.

,

Pringle

J. E.

1974

,

MNRAS

,

168

,

603

Crossref

Search ADS

Macías

E.

,

Guerra-Alvarado

O.

,

Carrasco-González

C.

,

Ribas

Á.

,

Espaillat

C. C.

,

Huang

J.

,

Andrews

S. M.

2021

,

A&A

,

648

,

A33

Crossref

Search ADS

Miyake

K.

,

Nakagawa

Y.

1993

,

Icarus

,

106

,

20

Crossref

Search ADS

Mori

S.

,

Bai

X.-N.

,

Okuzumi

S.

2019

,

ApJ

,

872

,

98

Crossref

Search ADS

Mori

S.

,

Okuzumi

S.

,

Kunitomo

M.

,

Bai

X.-N.

2021

,

ApJ

,

916

,

72

Crossref

Search ADS

Musiolik

G.

,

Wurm

G.

2019

,

ApJ

,

873

,

58

Crossref

Search ADS

Nakamoto

T.

,

Nakagawa

Y.

1994

,

ApJ

,

421

,

640

Crossref

Search ADS

Öberg

K. I.

,

Murray-Clay

R.

,

Bergin

E. A.

2011

,

ApJ

,

743

,

L16

Crossref

Search ADS

Ohno

K.

,

Ueda

T.

2021

,

A&A

,

651

,

L2

Crossref

Search ADS

Oka

A.

,

Nakamoto

T.

,

Ida

S.

2011

,

ApJ

,

738

,

141

Crossref

Search ADS

Okuzumi

S.

,

Momose

M.

,

Sirono

S.

,

Kobayashi

H.

,

Tanaka

H.

2016

,

ApJ

,

821

,

82

Crossref

Search ADS

Owen

J. E.

2020

,

MNRAS

,

495

,

3160

Crossref

Search ADS

Owen

J. E.

,

Armitage

P. J.

2014

,

MNRAS

,

445

,

2800

Crossref

Search ADS

Paardekooper

S. J.

,

Baruteau

C.

,

Crida

A.

,

Kley

W.

2010

,

MNRAS

,

401

,

1950

Crossref

Search ADS

Paardekooper

S. J.

,

Baruteau

C.

,

Kley

W.

2011

,

MNRAS

,

410

,

293

Crossref

Search ADS

Pavlyuchenkov

Y. N.

,

Maksimova

L. A.

,

Akimkin

V. V.

2022

,

Astron. Rep.

,

66

,

321

Crossref

Search ADS

Pinilla

P.

,

Pohl

A.

,

Stammler

S. M.

,

Birnstiel

T.

2017

,

ApJ

,

845

,

68

Crossref

Search ADS

Ros

K.

,

Johansen

A.

2013

,

A&A

,

552

,

A137

Crossref

Search ADS

Rosenfeld

K. A.

,

Andrews

S. M.

,

Hughes

A. M.

,

Wilner

D. J.

,

Qi

C.

2013

,

ApJ

,

774

,

16

Crossref

Search ADS

Ruden

S. P.

,

Pollack

J. B.

1991

,

ApJ

,

375

,

740

Crossref

Search ADS

Rybicki

G. B.

,

Lightman

A. P.

1979

,

Radiative Processes in Astrophysics

(

New York

:

Wiley

)

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Saito

E.

,

Sirono

S.

2011

,

ApJ

,

728

,

20

Crossref

Search ADS

Schoonenberg

D.

,

Ormel

C. W.

2017

,

A&A

,

602

,

A21

Crossref

Search ADS

Siebenmorgen

R.

,

Heymann

F.

2012

,

A&A

,

539

,

A20

Crossref

Search ADS

Sierra

A.

,

Lizano

S.

2020

,

ApJ

,

892

,

136

Crossref

Search ADS

Sirono

S.

2011

,

ApJ

,

733

,

L41

Crossref

Search ADS

Sirono

S.

,

Ueno

H.

2017

,

ApJ

,

841

,

36

Crossref

Search ADS

Stevenson

D. J.

,

Lunine

J. I.

1988

,

Icarus

,

75

,

146

Crossref

Search ADS

Turner

N. J.

,

Choukroun

M.

,

Castillo-Rogez

J.

,

Bryden

G.

2012

,

ApJ

,

748

,

92

Crossref

Search ADS

Ueda

T.

,

Flock

M.

,

Birnstiel

T.

2021

,

ApJ

,

914

,

L38

Crossref

Search ADS

Ueda

T.

,

Okuzumi

S.

,

Flock

M.

2017

,

ApJ

,

843

,

49

Crossref

Search ADS

Wada

K.

,

Tanaka

H.

,

Suyama

T.

,

Kimura

H.

,

Yamamoto

T.

2009

,

ApJ

,

702

,

1490

Crossref

Search ADS

Watanabe

S.

,

Lin

D. N. C.

2008

,

ApJ

,

672

,

1183

Crossref

Search ADS

Watanabe

S.

,

Nakagawa

Y.

,

Nakazawa

K.

1990

,

ApJ

,

358

,

282

Crossref

Search ADS

Weaver

E.

,

Isella

A.

,

Boehler

Y.

2018

,

ApJ

,

853

,

113

Crossref

Search ADS

Weidenschilling

S. J.

1977

,

MNRAS

,

180

,

57

Crossref

Search ADS

Whipple

F. L.

1972

, in

From Plasma to Planet

, ed.

Elvius

A.

(

New York

:

Wiley

),

211

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Wu

Y.

,

Lithwick

Y.

2021

,

ApJ

,

923

,

123

Crossref

Search ADS

Yang

C.-C.

,

Menou

K.

2010

,

MNRAS

,

402

,

2436

Crossref

Search ADS

Zhang

K.

,

Bergin

E. A.

,

Blake

G. A.

,

Cleeves

L. I.

,

Schwarz

K. R.

2017

,

Nature Astron.

,

1

,

0130

Crossref

Search ADS

Zhang

K.

,

Blake

G. A.

,

Bergin

E. A.

2015

,

ApJ

,

806

,

L7

Crossref

Search ADS

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
June 2022	9
July 2022	6
August 2022	13
September 2022	9
October 2022	6
November 2022	3
December 2022	9
January 2023	8
February 2023	5
March 2023	9
April 2023	10
May 2023	13
June 2023	17
July 2023	13
August 2023	20
September 2023	17
October 2023	21
November 2023	13
December 2023	15
January 2024	34
February 2024	19
March 2024	18
April 2024	17
May 2024	19
June 2024	19
July 2024	33
August 2024	20
September 2024	22
October 2024	15
November 2024	11
December 2024	13
January 2025	13
February 2025	15
March 2025	22
April 2025	5
May 2025	3

Article Contents

A global two-layer radiative transfer model for axisymmetric, shadowed protoplanetary disks

Abstract

1 Introduction

2 Local two-layer model

3 Global two-layer model

3.1 The irradiation surface

3.2 The thermal photosphere

3.3 Flux of the downward reprocessed starlight

3.4 Disk vertical structure

3.5 Disk temperature evolution

3.6 Numerical implementation

4 Validation with gapped disks

4.1 Model

4.2 Results

4.2.1 Mp = 0 model

4.2.2 Mp = 70M⊕ model

4.2.3 Mp = 200M⊕ model

4.2.4 The role of radial radiation transfer

5 Application to the thermal wave instability

5.1 Model

5.2 Results

5.3 Implications for dust evolution and planetesimal formation

6 Summary

Acknowledgements

Appendix 1. Fraction of starlight reprocessed downward

Appendix 2. Reprocessed starlight flux for plane-parallel disks

Appendix 3. Numerical treatment of the computational boundaries

Appendix 4. Emissivity of an isothermal slab

Appendix 5. The thermal wave instability: radial resolution dependence

Footnotes

References

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

4.2.1 M_p = 0 model

4.2.2 M_p = 70M_⊕ model

4.2.3 M_p = 200M_⊕ model