Analytic Expressions for the Inner-rim Structure of Passively Heated Protoplanetary Disks

Takahiro Ueda; Satoshi Okuzumi; Mario Flock

doi:10.3847/1538-4357/aa74b5

1. Introduction

The inner region of protoplanetary disks is the birthplace of rocky planetesimals and planets. One preferential site of rocky planetesimal formation is the inner edge of the so-called dead zone (e.g., Kretke et al. 2009). The dead zone is the location where magnetorotational instability (MRI; Balbus & Hawley 1998) is suppressed because of poor gas ionization (Gammie 1996). In the inner region of protoplanetary disks, the dead zone is likely to have an inner edge where the gas temperature T reaches ∼1000 K, above which the thermal ionization of the gas is effective enough to activate MRI (Gammie 1996; Desch & Turner 2015). Across the dead zone inner edge, the turbulent viscosity arising from MRI steeply decreases from inside out, resulting in a local maximum in the radial profile of the gas pressure (e.g., Dzyurkevich et al. 2010; Flock et al. 2016, 2017). The pressure maximum traps solid particles (Whipple 1972; Adachi et al. 1976) and may thereby facilitate dust growth (Brauer et al. 2008; Testi et al. 2014) and in situ planet formation (Tan et al. 2015). The dead zone edge may even trap planets by halting their inward migration (Masset et al. 2006). Because the location of the dead zone inner edge is determined by the gas temperature, understanding the temperature structure is important for understanding rocky planetesimal formation as well as the orbital architecture of inner planets.

The temperature structure of the innermost regions of protoplanetary disks has been extensively studied in the context of their observational appearance (see Dullemond & Monnier 2010 for review, and also Dullemond et al. 2001; Natta et al. 2001; Isella & Natta 2005). Recently, Flock et al. (2016) performed radiation hydrodynamical calculations of the the inner-rim structures of disks around Herbig Ae stars, taking into account the effects of starlight and viscous heating, evaporation and condensation of silicate grains, and gas opacity. They found that the inner disk consists of four distinct zones with different temperature profiles. According to their results, the dead zone inner edge occurs at the border between the two outermost regions, where the temperature drops steeply across 1000 K. Thus, this complex temperature structure potentially influences where rocky planetesimals and planets preferentially form. However, there has been no simple analytic model that reproduces the entire structure of the temperature profile. In addition, because Flock et al. (2016) performed simulations only for the disks around Herbig Ae stars, it has been unclear how the inner-rim structures depend on stellar parameters.

In this paper, we analytically derive the temperature and dust-to-gas mass ratio profile for the inner region of disks and reveal how these structures depend on the central star parameters based on the results from Flock et al. (2016). Our analysis and the analytic solutions are shown in Section 2. We discuss the implications for the planetesimal and planet formation and the disk observation in Section 3. A summary is in Section 4.

2. Analytic Solutions for the Inner-rim Structure

We derive analytic expressions for the inner-rim structure of protoplanetary disks based on the recent numerical simulations conducted by Flock et al. (2016). They found that the inner region of protoplanetary disks can be divided into four regions (regions A, B, C, and D) in terms of the temperature profile. Figure 1 schematically shows the inner region of protoplanetary disks. In the following subsections, we describe how structures such as the temperature and the dust-to-gas mass ratio are determined in each region.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Schematic illustration of the inner structure of a passively irradiated protoplanetary disk. The disk received the starlight on the surface where the visual optical depth ${\tau }_{* }$ measured along the starlight reaches unity (solid curve), and receives the reprocessed light from the ${\tau }_{* }=1$ surface on the deeper surface where the infrared optical depth ${\tau }_{\mathrm{IR}}$ measured vertically (z-direction) from above reaches unity (dashed curve). The inner region can be divided into four regions (regions A, B, C, and D). Region A: an optically thin dust-free region; region B: an optically thin dust halo region where the temperature is nearly equal to the evaporation temperature; region C: an optically thick dust condensation front irradiated by the central star, with a high incident angle; region D: an optically thick disk illuminated by not only the central star but also the dust condensation front.
Download figure:
Standard image High-resolution image

**Figure 1.** Schematic illustration of the inner structure of a passively irradiated protoplanetary disk. The disk received the starlight on the surface where the visual optical depth ${\tau }_{* }$ measured along the starlight reaches unity (solid curve), and receives the reprocessed light from the ${\tau }_{* }=1$ surface on the deeper surface where the infrared optical depth ${\tau }_{\mathrm{IR}}$ measured vertically (z-direction) from above reaches unity (dashed curve). The inner region can be divided into four regions (regions A, B, C, and D). Region A: an optically thin dust-free region; region B: an optically thin dust halo region where the temperature is nearly equal to the evaporation temperature; region C: an optically thick dust condensation front irradiated by the central star, with a high incident angle; region D: an optically thick disk illuminated by not only the central star but also the dust condensation front.
Download figure:
Standard image High-resolution image

2.1. Region A: An Optically Thin Dust-free Region

Region A is the innermost region where the temperature is above the evaporation temperature of dust. This region is therefore free of dust and is optically thin to the starlight. In an optically thin region, the temperature profile is generally given by

$\begin{eqnarray}&&T={\epsilon }^{-1/4}{\left(\displaystyle \frac{{R}_{* }}{2r}\right)}^{1/2}{T}_{* },\end{eqnarray} \tag{ 1 }$

where r is the distance from the central star, and R_* and T_* are the stellar radius and temperature, respectively. The dimensionless quantity is the ratio between the emission and absorption efficiencies of the disk gas, including the contribution from dust,

$\begin{eqnarray}&&\epsilon \equiv \displaystyle \frac{{\kappa }_{{\rm{g}}}+{f}_{{\rm{d}}2{\rm{g}}}{\kappa }_{{\rm{d}}}(T)}{{\kappa }_{{\rm{g}}}+{f}_{{\rm{d}}2{\rm{g}}}{\kappa }_{{\rm{d}}}({T}_{* })},\end{eqnarray} \tag{ 2 }$

where ${f}_{{\rm{d}}2{\rm{g}}}$ is the dust-to-gas mass ratio, ${\kappa }_{{\rm{g}}}$ is the Planck mean gas opacity (here assumed to be independent of the temperature), and ${\kappa }_{{\rm{d}}}({T}_{* })$ and ${\kappa }_{{\rm{d}}}(T)$ are the Planck mean dust opacities at the stellar and disk temperatures, respectively. Flock et al. (2016) adopted ${\kappa }_{{\rm{g}}}={10}^{-4}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},{\kappa }_{{\rm{d}}}({T}_{* })=2100\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ , and ${\kappa }_{{\rm{d}}}(T)=700\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ .

For region A, $\epsilon \approx 1$ because ${f}_{{\rm{d}}2{\rm{g}}}\approx 0$ . Therefore, the temperature profile in region A, ${T}_{{\rm{A}}}$ , is given by

$\begin{eqnarray}&&{T}_{{\rm{A}}}={\left(\displaystyle \frac{{R}_{* }}{2r}\right)}^{1/2}{T}_{* }.\end{eqnarray} \tag{ 3 }$

Flock et al. (2016) showed that Equation (3) accurately reproduces the temperature profile in region A (see the top panel of their Figure 1). The real gas temperature might be more difficult to understand because it depends on gas opacity (Dullemond & Monnier 2010; Hirose 2015), which is assumed to be independent of wavelength in this paper.

2.2. Region B: An Optically Thin Dust Halo Region

Region B is the location where dust starts to condense but is still optically thin. This region, which Flock et al. (2016) called the dust halo, has a uniform temperature that is nearly equal to the evaporation temperature ${T}_{\mathrm{ev}}$ (see Figure 1 of Flock et al. 2016). These features can be understood by noting that the condensing dust acts as a thermostat: if T falls below ${T}_{\mathrm{ev}}$ , dust starts to condense, but the condensed dust pushes the temperature back up because the emission-to-absorption ratio of the dust is lower than that of the gas. Thus, dust evaporation and condensation regulate the temperature in this region to $\sim {T}_{\mathrm{ev}}$ . Strictly speaking, the temperature of region B obtained by Flock et al. (2016) is higher than ${T}_{\mathrm{ev}}$ by about 100 K. However, this is due to their numerical treatment of the dust evaporation, where the dust-to-gas ratio is given by a smoothed step function of T with a smoothing width of 100 K.

$T\sim {T}_{\mathrm{ev}}$ in region B can be used to predict the spatial distribution of the dust-to-gas mass ratio ${f}_{{\rm{d}}2{\rm{g}}}$ in this region. Substituting Equation (2) into Equation (1) and solving the equation with respect to ${f}_{{\rm{d}}2{\rm{g}}}$ , we obtain

$\begin{eqnarray}&&{f}_{{\rm{d}}2{\rm{g}}}=\displaystyle \frac{{\kappa }_{{\rm{g}}}(4{r}^{2}{T}_{{\rm{B}}}^{4}-{R}_{* }^{2}{T}_{* }^{4})}{{\kappa }_{{\rm{d}}}({T}_{* }){R}_{* }^{2}{T}_{* }^{4}-4{\kappa }_{{\rm{d}}}({T}_{{\rm{B}}}){r}^{2}{T}_{{\rm{B}}}^{4}},\end{eqnarray} \tag{ 4 }$

where ${T}_{{\rm{B}}}$ stands for the temperature in region B. Assuming that ${T}_{{\rm{B}}}$ is constant, Equation (4) determines ${f}_{{\rm{d}}2{\rm{g}}}$ as a function of r. Figure 2 compares Equation (4) with the radial profile of ${f}_{{\rm{d}}2{\rm{g}}}$ at the midplane taken from the radiation hydrostatic disk model S100 of Flock et al. (2016; see the bottom panel of their Figure 1. The stellar temperature, radius, mass, and luminosity are set to ${T}_{* }={\rm{10,000}}\,{\rm{K}}$ , ${R}_{* }=2.5{R}_{\odot }=0.0116\,\mathrm{au}$ , ${M}_{* }=2.5{M}_{\odot }$ and ${L}_{* }=56{L}_{\odot }$ , respectively. See also their Table 1 for details). Equation (4) perfectly reproduces the results of Flock et al. (2016) when ${T}_{{\rm{B}}}$ is set to ${T}_{\mathrm{ev}}+100\,{\rm{K}}\approx 1470\,{\rm{K}}$ , where we have used that ${T}_{\mathrm{ev}}$ in region B is $\approx 1370\,{\rm{K}}$ in their calculation.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Dust-to-gas mass ratio profile ${f}_{{\rm{d}}2{\rm{g}}}$ at the midplane in region B. The red and black solid lines show our analytic solution (Equation (4)) and the radiation hydrostatic disk model S100 of Flock et al. (2016) with the highest resolution, respectively.
Download figure:
Standard image High-resolution image

As we can see from Figure 2, the inner and outer edges of region B correspond to the locations where ${f}_{{\rm{d}}2{\rm{g}}}$ given by Equation (4) goes to zero and infinity, respectively. The radii of these locations are given by

$\begin{eqnarray}{R}_{\mathrm{AB}} & = & \displaystyle \frac{1}{2}{\left(\displaystyle \frac{{T}_{* }}{{T}_{{\rm{B}}}}\right)}^{2}{R}_{* }\\ & = & 0.27{\left(\displaystyle \frac{{T}_{* }}{{10}^{4}{\rm{K}}}\right)}^{2}{\left(\displaystyle \frac{{T}_{{\rm{B}}}}{1470{\rm{K}}}\right)}^{-2}\left(\displaystyle \frac{{R}_{* }}{2.5{R}_{\odot }}\right)\,\mathrm{au}\end{eqnarray} \tag{ 5 }$

and

$\begin{eqnarray}{R}_{\mathrm{BC}} & = & \displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\kappa }_{{\rm{d}}}({T}_{* })}{{\kappa }_{{\rm{d}}}({T}_{{\rm{B}}})}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{* }}{{T}_{{\rm{B}}}}\right)}^{2}{R}_{* }\\ & = & 0.46{\left(\displaystyle \frac{{\kappa }_{{\rm{d}}}({T}_{* })/{\kappa }_{{\rm{d}}}({T}_{{\rm{B}}})}{3}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{* }}{{10}^{4}{\rm{K}}}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{{T}_{{\rm{B}}}}{1470{\rm{K}}}\right)}^{-2}\left(\displaystyle \frac{{R}_{* }}{2.5{R}_{\odot }}\right)\,\mathrm{au},\end{eqnarray} \tag{ 6 }$

respectively. We have expressed the inner radius of region B as ${R}_{\mathrm{AB}}$ because it also stands for the outer radius of the dust-free zone (region A). Similarly, ${R}_{\mathrm{BC}}$ also stands for the inner radius of the optically thick inner-rim (region C). We note that our expression for ${R}_{\mathrm{BC}}$ is equivalent to that of Monnier & Millan-Gabet (2002).

Equations (5) and (6) give the orbital radii of the boundaries of region B if ${T}_{{\rm{B}}}(\approx \,{T}_{\mathrm{ev}})$ is given. In fact, ${T}_{\mathrm{ev}}$ weakly depends on the gas density, and hence on r. Therefore, if one wants to determine ${R}_{\mathrm{AB}}$ and ${R}_{\mathrm{BC}}$ precisely, one has to solve Equations (5) and (6) simultaneously with the equation for ${T}_{\mathrm{ev}}$ as a function of r. This is demonstrated in the Appendix.

As we showed, a dust halo naturally appears. The dust halo is a main contributor of the near-infrared emission in the spectral energy distribution of disks around Herbig Ae/Be stars (e.g., Vinković et al. 2006).

2.3. Region C: An Optically Thick Dust Condensation Front

At $r\approx {R}_{\mathrm{BC}}$ , dust grains near the midplane fully condense (see Figure 2), and the visual optical depth ${\tau }_{* }$ measured along the line directly from the star exceeds unity at the midplane. Further out, stellar irradiation is absorbed at a height where ${\tau }_{* }=1$ , and the midplane temperature is determined by the reprocessed, infrared radiation from dust grains lying at the ${\tau }_{* }=1$ surface. Figure 3 schematically illustrates the surface structure of regions C and D. As we show below, region C can be regarded as the zone where the ${\tau }_{* }=1$ surface coincides with the condensation front; the surface structure is approximately determined by the condition that the surface temperature ${T}_{{\rm{s}},{\rm{C}}}$ must be equal to the characterisctic temperature $\sim 1200\,{\rm{K}}$ .

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Schematic illustration of the geometric and radiative structures of regions C and D. The stellar radiation F_* determines the temperature of the optically thin layer above the ${\tau }_{* }=1$ surface. The dust on the ${\tau }_{* }=1$ surface reprocesses the starlight into infrared radiation ${F}_{\mathrm{IR}}$ , and approximately half of the dust emission gets absorbed by the midplane region of the disk. In region C, dust condensation/evaporation controls the location of the ${\tau }_{* }=1$ surface.
Download figure:
Standard image High-resolution image

**Figure 3.** Schematic illustration of the geometric and radiative structures of regions C and D. The stellar radiation F_* determines the temperature of the optically thin layer above the ${\tau }_{* }=1$ surface. The dust on the ${\tau }_{* }=1$ surface reprocesses the starlight into infrared radiation ${F}_{\mathrm{IR}}$ , and approximately half of the dust emission gets absorbed by the midplane region of the disk. In region C, dust condensation/evaporation controls the location of the ${\tau }_{* }=1$ surface.
Download figure:
Standard image High-resolution image

The radial extent of region C is of particular interest because the outer edge of this region sets the inner edge of the dead zone (Flock et al. 2016, 2017). We will see in Section 2.4 that the border of regions C and D is related to how the height of the ${\tau }_{* }=1$ surface, z_*, depends on disk radius. Here, we try to estimate this height by considering the energy balance between stellar irradiation on the ${\tau }_{* }=1$ surface and the thermal emission from the surface,

$\begin{eqnarray}&&\displaystyle \frac{{L}_{* }}{4\pi {R}^{2}}\sin \theta ={\left(\displaystyle \frac{{C}_{\mathrm{bw}}}{4}\right)}^{-1}\epsilon {\sigma }_{\mathrm{SB}}{T}_{{\rm{s}},{\rm{C}}}^{4},\end{eqnarray} \tag{ 7 }$

where R is the radial distance on the midplane, ${\sigma }_{\mathrm{SB}}$ is the Stefan–Boltzmann constant, and θ is the angle between the starlight and the disk surface (the so-called grazing angle). ${C}_{\mathrm{bw}}$ is a backwarming factor (e.g., Dullemond et al. 2001), which is defined as the ratio of a full $4\pi$ solid angle to the solid angle subtended by the empty sky seen by particles on the surface. The backwarming is heating by surrounding dust particles. If a dust particle radiates into the sky area covered by nearby dust particles, a similar amount of energy is returned from surrounding particles (Kama et al. 2009). We assume ${C}_{\mathrm{bw}}=4$ , which means that the ${\tau }_{* }=1$ surface can be approximated as a flat wall. In region C, $\epsilon =1/3$ because dust totally condenses. The grazing angle θ is related to the surface height ${z}_{* }(R)$ as (Chiang & Goldreich 1997; Tanaka et al. 2005)

$\begin{eqnarray}\theta & = & \arcsin \left(\displaystyle \frac{4{R}_{* }}{3\pi R}\right)+\arctan \left(\displaystyle \frac{{z}_{* }}{R}\displaystyle \frac{d\mathrm{ln}{z}_{* }}{d\mathrm{ln}R}\right)-\arctan \left(\displaystyle \frac{{z}_{* }}{R}\right)\\ & \approx & \displaystyle \frac{4{R}_{* }}{3\pi R}+\displaystyle \frac{{z}_{* }}{R}\displaystyle \frac{d\mathrm{ln}{z}_{* }}{d\mathrm{ln}R}-\displaystyle \frac{{z}_{* }}{R},\end{eqnarray} \tag{ 8 }$

where the second expression assumes $\theta \ll 1$ . As mentioned earlier, the ${\tau }_{* }=1$ surface in region C coincides with the condensation front, and the surface temperature ${T}_{{\rm{s}},{\rm{C}}}$ is approximately equal to 1200 K according to the simulations by Flock et al. (2016; see their Figure 2). Substituting Equation (8) and ${T}_{{\rm{s}},{\rm{C}}}\approx \mathrm{constant}$ into Equation (7) and solving the differential equation with respect to z_*, we obtain the surface height in region C as

$\begin{eqnarray}&&{z}_{* ,{\rm{C}}}=\displaystyle \frac{1}{6{R}_{* }^{2}}{\left(\displaystyle \frac{{T}_{{\rm{s}},{\rm{C}}}}{{T}_{* }}\right)}^{4}({R}^{3}-{{RR}}_{0}^{2})+\displaystyle \frac{4{R}_{* }}{3\pi {R}_{0}}(R-{R}_{0}),\end{eqnarray} \tag{ 9 }$

where R₀ is the radius at which ${z}_{* ,{\rm{C}}}=0$ . Figure 4 shows the location of the ${\tau }_{* }=1$ plane obtained by Equation (9) and that from 3D simulations by Flock et al. (2016). In Figure 4, we set ${R}_{0}=0.4\,\mathrm{au}$ , which explains Flock et al. (2016) better than ${R}_{0}={R}_{\mathrm{BC}}\approx 0.46\,\mathrm{au}$ . Our analytic solution (Equation (9)) reproduces the shape of the condensation front obtained by Flock et al. (2016).

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Locations of the irradiation ( ${\tau }_{* }=1$ ) surface for the disk model S100 of Flock et al. (2016) obtained from our analytic expressions (solid curves) and from the radiation hydrostatic calculation by Flock et al. (2016; dashed curve). The red solid curve shows the location of the irradiation surface in region C from our Equation (9) with ${R}_{0}=0.4\,\mathrm{au}$ . The blue and green curves show ${z}_{* ,{\rm{D}}}=3.6{h}_{{\rm{g}}}$ and $4.8{h}_{{\rm{g}}}$ , which well approximate the height of the irradiation surface in the inner and outer parts of region D, respectively. The curves are overlaid on the 2D temperature map obtained from the radiation hydrostatic calculation of Flock et al. (2016).
Download figure:
Standard image High-resolution image

In order to obtain the midplane temperature in region C, ${T}_{\mathrm{mid},{\rm{C}}}$ , we assume that half of the infrared emission from the surface layer comes into the interior. Then we have

$\begin{eqnarray}&&{T}_{\mathrm{mid},{\rm{C}}}={2}^{-1/4}{T}_{{\rm{s}},{\rm{C}}}\approx 1009\,{\rm{K}}.\end{eqnarray} \tag{ 10 }$

2.4. Region D: An Optically Thick Region

Region D is the outermost, coldest region where dust condenses at all heights. Therefore, the temperature structure of this region is essentially the same as that of the classical two-layer disk model by Chiang & Goldreich (1997) and Kusaka et al. (1970). Specifically, the surface temperature profile in this region follows Equations (7) with (8), $\epsilon =1/3$ , and ${C}_{\mathrm{bw}}=4$ , i.e.,

$\begin{eqnarray}&&{T}_{{\rm{s}},{\rm{D}}}^{4}=3{\left(\displaystyle \frac{R}{{R}_{* }}\right)}^{-2}\left[\displaystyle \frac{4{R}_{* }}{3\pi R}+\displaystyle \frac{{z}_{* ,{\rm{D}}}}{R}\left(\displaystyle \frac{d\mathrm{ln}{z}_{* ,{\rm{D}}}}{d\mathrm{ln}R}-1\right)\right]{T}_{* }^{4},\end{eqnarray} \tag{ 11 }$

where ${z}_{* ,{\rm{D}}}$ is the surface height in region D. Assuming that ${z}_{* ,{\rm{D}}}$ is proportional to the gas scale height ${h}_{{\rm{g}}}={c}_{{\rm{s}}}/{{\rm{\Omega }}}_{{\rm{K}}}$ , where ${c}_{{\rm{s}}}$ is the sound speed and ${{\rm{\Omega }}}_{{\rm{K}}}$ is Keplerian angular velocity, we obtain (Kusaka et al. 1970)

$\begin{eqnarray}&&{T}_{{\rm{s}},{\rm{D}}}^{4}={T}_{1}^{4}+{T}_{2}^{4},\end{eqnarray} \tag{ 12 }$

where

$\begin{eqnarray}{T}_{1} & = & {\left(\displaystyle \frac{4}{\pi }\right)}^{1/4}{\left(\displaystyle \frac{R}{{R}_{* }}\right)}^{-3/4}{T}_{* }\\ & = & 333{\left(\displaystyle \frac{R}{1\mathrm{au}}\right)}^{-3/4}{\left(\displaystyle \frac{{R}_{* }}{2.5{{R}}_{\odot }}\right)}^{3/4}\left(\displaystyle \frac{{T}_{* }}{{10}^{4}\,{\rm{K}}}\right)\,{\rm{K}}\end{eqnarray} \tag{ 13 }$

and

$\begin{eqnarray}{T}_{2} & = & {\left(\displaystyle \frac{6}{7}\right)}^{2/7}{\left(\displaystyle \frac{{z}_{* ,{\rm{D}}}}{{h}_{{\rm{g}}}}\right)}^{2/7}{\left(\displaystyle \frac{R}{{R}_{* }}\right)}^{-3/7}{\left(\displaystyle \frac{{k}_{{\rm{B}}}{T}_{* }{R}_{* }}{{{mGM}}_{* }}\right)}^{1/7}{T}_{* }\\ & = & 413{\left(\displaystyle \frac{{z}_{* ,{\rm{D}}}}{{h}_{{\rm{g}}}}\right)}^{2/7}{\left(\displaystyle \frac{R}{1\mathrm{au}}\right)}^{-3/7}{\left(\displaystyle \frac{{R}_{* }}{2.5{{R}}_{\odot }}\right)}^{4/7}\\ & & \times {\left(\displaystyle \frac{{M}_{* }}{2.5{M}_{\odot }}\right)}^{-1/7}{\left(\displaystyle \frac{{T}_{* }}{{10}^{4}{\rm{K}}}\right)}^{8/7}\,{\rm{K}},\end{eqnarray} \tag{ 14 }$

where ${k}_{{\rm{B}}}$ is the Boltzmann constant, G is the gravitational constant, and $m=4\times {10}^{-24}\,{\rm{g}}$ is the mean molecular mass of the disk gas. As in region C, the midplane temperature of region D is given by

$\begin{eqnarray}&&{T}_{\mathrm{mid},{\rm{D}}}={2}^{-1/4}{T}_{{\rm{s}},{\rm{D}}}.\end{eqnarray} \tag{ 15 }$

The border between regions C and D can be defined as the position where the ${\tau }_{* }=1$ surfaces in the two regions intersect. Equating ${z}_{* ,{\rm{C}}}$ given by Equation (9) with ${z}_{* ,{\rm{D}}}=({z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}){h}_{{\rm{g}}}=({z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}})\sqrt{{k}_{{\rm{B}}}{T}_{\mathrm{mid},{\rm{D}}}{R}^{3}/{{mGM}}_{* }}$ , we obtain the transcendental equation for the radius ${R}_{\mathrm{CD}}$ of the border,

$\begin{eqnarray}&&\displaystyle \frac{1}{6{R}_{* }^{2}}{\left(\displaystyle \frac{{T}_{{\rm{s}},{\rm{C}}}}{{T}_{* }}\right)}^{4}({R}_{\mathrm{CD}}^{2}-{R}_{\mathrm{BC}}^{2})+\displaystyle \frac{4{R}_{* }}{3\pi {R}_{\mathrm{BC}}{R}_{\mathrm{CD}}}({R}_{\mathrm{CD}}-{R}_{\mathrm{BC}})\\ &&\quad =\,\left(\displaystyle \frac{{z}_{* ,{\rm{D}}}}{{h}_{{\rm{g}}}}\right)\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{mid},\mathrm{CD}}{R}_{\mathrm{CD}}}{{{mGM}}_{* }}},\end{eqnarray} \tag{ 16 }$

where the subscript CD stands for the value at $R={R}_{\mathrm{CD}}$ and we have used ${R}_{0}\approx {R}_{\mathrm{BC}}$ (see Section 2.3). In order to derive a closed-form expression for ${R}_{\mathrm{CD}}$ , we simplify Equation (16) as follows. First, we neglect the second term on the left side of Equation (16) because for T-Tauri and Herbig stars the inner-rim radius is much larger than the stellar radius. Second, we replace the factor ${T}_{\mathrm{mid},\mathrm{CD}}^{1/2}={({2}^{-1}{T}_{{\rm{s}},\mathrm{CD}}^{4})}^{1/8}$ on the right side of Equation (16) with ${({T}_{2,\mathrm{CD}}^{4})}^{1/8}$ , where ${T}_{2,\mathrm{CD}}$ is the value of T₂ at $R={R}_{\mathrm{CD}}$ , because ${T}_{1}\sim {T}_{2}$ at R ∼ 0.1–1 au (see Equations (13) and (14)). Third, we neglect the ${R}_{\mathrm{CD}}$ -dependence of the right side, by substituting ${R}_{\mathrm{CD}}=2{R}_{\mathrm{BC}}$ on the right side only, because the dependence is generally weak ( $\sqrt{{R}_{\mathrm{CD}}{T}_{\mathrm{mid},\mathrm{CD}}}\,\propto {R}_{\mathrm{CD}}^{1/8}$ for ${T}_{\mathrm{mid},\mathrm{CD}}\approx {2}^{-1/4}{T}_{1,\mathrm{CD}}$ , and $\sqrt{{R}_{\mathrm{CD}}{T}_{\mathrm{mid},\mathrm{CD}}}\propto {R}_{\mathrm{CD}}^{2/7}$ for ${T}_{\mathrm{mid},\mathrm{CD}}\approx {2}^{-1/4}{T}_{2,\mathrm{CD}}$ ). With these simplifications, Equation (16) can be analytically solved with respect to ${R}_{\mathrm{CD}}$ , and we thus obtain a closed-form expression for ${R}_{\mathrm{CD}}$ ,

$\begin{eqnarray}&&{R}_{\mathrm{CD}}=\sqrt{1+{\rm{\Gamma }}}{R}_{\mathrm{BC}},\end{eqnarray} \tag{ 17 }$

where

$\begin{eqnarray}{\rm{\Gamma }} & \equiv & 3.1{\left(\displaystyle \frac{{R}_{\mathrm{BC}}}{0.46\mathrm{au}}\right)}^{-12/7}{\left(\displaystyle \frac{{z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}}{4.8}\right)}^{8/7}{\left(\displaystyle \frac{{T}_{* }}{{10}^{4}{\rm{K}}}\right)}^{32/7}\\ & & \times {\left(\displaystyle \frac{{M}_{* }}{2.5{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{* }}{2.5{R}_{\odot }}\right)}^{16/7}.\end{eqnarray} \tag{ 18 }$

If we use Equation (6) to eliminate ${R}_{\mathrm{BC}}$ from Γ, we obtain

$\begin{eqnarray}&&{\rm{\Gamma }}=3.1{\left(\displaystyle \frac{{z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}}{4.8}\right)}^{8/7}{\left(\displaystyle \frac{{L}_{* }}{56{L}_{\odot }}\right)}^{2/7}{\left(\displaystyle \frac{{M}_{* }}{2.5{M}_{\odot }}\right)}^{-1/2},\end{eqnarray} \tag{ 19 }$

where we set ${\kappa }_{{\rm{d}}}({T}_{* })/{\kappa }_{{\rm{d}}}({T}_{{\rm{B}}})=1/3$ and ${T}_{{\rm{B}}}=1470\,{\rm{K}}$ . The resultant equations (Equations (17) and (19)) justify the validity of the approximation of ${R}_{\mathrm{CD}}=2{R}_{\mathrm{BC}}$ . Our simplifications introduce an error of less than $\,10$ % in the estimate of ${R}_{\mathrm{CD}};$ if we solve Equation (16) exactly, we obtain ${R}_{\mathrm{CD}}=0.85\,\mathrm{au}$ with ${T}_{* }={\rm{10,000}}\,{\rm{K}}$ , ${R}_{* }=2.5{R}_{\odot }$ , ${M}_{* }=2.5{M}_{\odot }$ , and ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=4.8$ , while Equation (17) leads to ${R}_{\mathrm{CD}}=0.92\,\mathrm{au}$ . As seen in Figure 4, the ratio ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}$ is a function of the distance from the central star: the ratio ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}$ is about 3.6 at the outer part and 4.8 at the inner part of region D, due to the extra heating by the hot rim surface discussed below. Therefore, when we estimate the location of ${R}_{\mathrm{CD}}$ , it is better to use the higher value for ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}$ (≈ 4.8). On the other hand, when we estimate the temperature in region D, it is better to use the lower value for ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}$ (≈ 3.6).

In Figure 5, we compare our analytic temperature profile with the profile derived from the radiative hydrostatic disk model S100 of Flock et al. (2016). We find that our analytic profile reproduces the result of Flock et al. (2016), except in the inner part of region D (R ≈ 0.8–1 au), where our Equation (15) underestimates the midplane temperature. This implies that some sources other than the central star heat the inner part of region D.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Radial profile of the midplane temperature ${T}_{\mathrm{mid}}$ for the disk model S100 of Flock et al. (2016). The black solid line is from the radiative hydrostatic calculation by Flock et al. (2016; their Figure 1), whereas the red solid line is from our analytic expressions: Equation (3) for region A, $T={T}_{\mathrm{eva}}=1470\,{\rm{K}}$ for region B, Equation (10) for region C, and Equation (15) with ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=3.6$ for region D. Boundaries are determined by Equations (5), (6), and (17). For reference, the temperature profiles for optically thin disks with $\epsilon =1$ and 1/3 are shown by the blue and purple dashed lines, respectively, and the profile for a smooth optically thick disk (Equation (15)) is shown by the green dashed line.
Download figure:
Standard image High-resolution image

As we show below, this discrepancy can be resolved by taking into account the irradiation by the hot rim surface located in region C as shown in Figure 1. For simplicity, we approximate the ${\tau }_{* }=1$ surface of region C by a flat surface truncated at $R={R}_{\mathrm{CD}}$ , and assume that the ${\tau }_{\mathrm{IR}}=1$ surface is located at a distance h below the ${\tau }_{* }=1$ surface (see Figure 6). At any radial position R on the ${\tau }_{\mathrm{IR}}$ surface, the balance between radiative cooling on the ${\tau }_{\mathrm{IR}}=1$ surface and irradiation heating by the hot rim surface can be written as

$\begin{eqnarray}&&{\sigma }_{\mathrm{SB}}{T}_{\mathrm{mid},\ {\rm{D}}}^{4}={\int }_{-\infty }^{{R}_{\mathrm{CD}}}\displaystyle \frac{{\sigma }_{\mathrm{SB}}{T}_{{\rm{s}},{\rm{C}}}^{4}}{2\pi l}\displaystyle \frac{h}{l}{dR}^{\prime} \end{eqnarray} \tag{ 20 }$

with $l=\sqrt{{h}^{2}+{(R-R^{\prime} )}^{2}}$ , where we have used the fact that in the absence of internal heat sources, the temperature at the ${\tau }_{\mathrm{IR}}=1$ is equal to the midplane surface. Performing the integration, we obtain

$\begin{eqnarray}&&{T}_{\mathrm{mid},{\rm{D}}}={\left\{\displaystyle \frac{1}{2\pi }\left[{\tan }^{-1}\left(\displaystyle \frac{{R}_{\mathrm{CD}}-R}{h}\right)+\displaystyle \frac{\pi }{2}\right]\right\}}^{1/4}{T}_{{\rm{s}},{\rm{C}}}.\end{eqnarray} \tag{ 21 }$

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Schematic illustration of the irradiation by the condensation front at the vicinity of the boundary between regions C and D. A line element with a width ${dR}^{\prime}$ at a position $R^{\prime}$ on the truncated ${\tau }_{* }=1$ surface irradiates in the infrared wavelength with temperature ${T}_{{\rm{s}},{\rm{C}}}$ . The emission travels along the distance l and gets absorbed by the position R on the ${\tau }_{\mathrm{IR}}=1$ surface located at distance h below the ${\tau }_{* }=1$ surface.
Download figure:
Standard image High-resolution image

${dR}^{\prime} $ — **Figure 6.** Schematic illustration of the irradiation by the condensation front at the vicinity of the boundary between regions C and D. A line element with a width ${dR}^{\prime}$ at a position $R^{\prime}$ on the truncated ${\tau }_{* }=1$ surface irradiates in the infrared wavelength with temperature ${T}_{{\rm{s}},{\rm{C}}}$ . The emission travels along the distance l and gets absorbed by the position R on the ${\tau }_{\mathrm{IR}}=1$ surface located at distance h below the ${\tau }_{* }=1$ surface.
Download figure:
Standard image High-resolution image

Strictly speaking, the ${\tau }_{* }=1$ surface of region C is not parallel to the ${\tau }_{\mathrm{IR}}=1$ surface of region D. From Equations (9) and (15), the inclinations of ${\tau }_{* }=1$ surfaces in regions C and D with respect to the midplane at $R={R}_{\mathrm{CD}}$ can be estimated as 30^◦ and 10^◦, respectively. Hence, if we assume that ${\tau }_{\mathrm{IR}}=1$ surfaces in regions C and D are parallel to the ${\tau }_{* }=1$ surfaces in regions C and D, respectively, the angle made by the ${\tau }_{* }=1$ surface in region C and the ${\tau }_{\mathrm{IR}}=1$ surface in region D is ≈20°. This deviation from the plane-parallel approximation introduces an error of ∼20% in temperature. Nevertheless, our analytic formula (Equation (21)) based on the approximation well reproduces the result of Flock et al. (2016) with an error much less than ∼20%. This might be because the actual position of ${\tau }_{\mathrm{IR}}=1$ surface in region D is not parallel to the ${\tau }_{* }=1$ surface in region D and should be determined by the optical depth along the ray emitted by each position on the ${\tau }_{* }=1$ surface in region C. If one needs a more accurate temperature structure, detailed radiative transfer simulations would be needed.

Figure 7 shows the analytic temperature profile refined with Equation (21). We set $h=0.05{R}_{\mathrm{CD}}$ , which is a bit larger than the typical gas scale height and ${R}_{\mathrm{CD}}=0.86\,\mathrm{au}$ , which explains Flock et al. (2016) better than the value estimated from Equation (17) (≈0.83 au with ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=3.6$ ). The larger value of h leads to the larger radial extent of the region where the additional heating by the hot rim is important (region ${\rm{D}}^{\prime}$ ). Here, we use Equations (21) or (15), whichever is greater, for the temperature in region D. We find that the refined model reproduces the temperature profile of Flock et al. (2016) in the vicinity of the boundary between regions C and D.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Same as Figure 5, except the analytic temperature profile (red line) now uses Equation (21), with $h=0.05{R}_{\mathrm{CD}}$ for the midplane temperature in the inner part of region D, region ${\rm{D}}^{\prime}$ . We set ${R}_{\mathrm{CD}}=0.86\,\mathrm{au}$ , which explains Flock et al. (2016) better than the value estimated from Equation (17) ( $0.83\,\mathrm{au}$ with ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=3.6$ ).
Download figure:
Standard image High-resolution image

3. Discussion

3.1. The Position of the Dead Zone Inner Edge

One of the goals of this work is to determine the position of the dead zone inner edge, ${R}_{\mathrm{DIB}}$ . Our results show that at the boundary between regions C and D, the temperature steeply decreases from ∼1000 K to ∼600 K. According to Desch & Turner (2015), the position of the dead zone inner edge should be located where $T\sim 1000\,{\rm{K}}$ . Therefore, the boundary between regions C and D is thought to be the dead zone inner edge. From Equation (17), ${R}_{\mathrm{DIB}}$ can be estimated as $0.09\,\mathrm{au}$ , with ${T}_{* }=5000\,{\rm{K}}$ and ${M}_{* }=1{M}_{\odot }$ , and $0.8\,\mathrm{au}$ with ${T}_{* }={\rm{10,000}}\,{\rm{K}}$ and ${M}_{* }=2.5{M}_{\odot }$ , with ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=4.8$ , ${\kappa }_{{\rm{d}}}({T}_{* })/{\kappa }_{{\rm{d}}}({T}_{{\rm{B}}})=1/3$ , and ${T}_{{\rm{B}}}=1470\,{\rm{K}}$ , while the classical optically thick temperature profile (Equation (15)) leads to ${R}_{\mathrm{DIB}}\approx 0.03\,\mathrm{au}$ , with ${T}_{* }=5000\,{\rm{K}}$ and ${M}_{* }=1{M}_{\odot }$ and ${R}_{\mathrm{DIB}}\approx 0.3\,\mathrm{au}$ with ${T}_{* }={\rm{10,000}}\,{\rm{K}}$ and ${M}_{* }=2.5{M}_{\odot }$ . Therefore, the location of the dead zone inner edge estimated from our model is 2–3 times farther out than that estimated from the classical optically thick temperature profile. This is because the temperature in region C is much higher than that estimated from the classical model due to the stellar irradiation on the condensation front with a high incident angle.

3.2. The Effect of Viscous Heating

We have focused on passively irradiated disks where the effect of viscous heating on the temperature structure is negligibly small. In the opposite case where viscous heating dominates over stellar irradiation, the midplane temperature ${T}_{\mathrm{mid},\mathrm{vis}}$ is given by

$\begin{eqnarray}&&{T}_{\mathrm{mid},\mathrm{vis}}={\left(1+\displaystyle \frac{3{\tau }_{\mathrm{IR}}}{4}\right)}^{1/4}{T}_{{\rm{s}},\mathrm{vis}},\end{eqnarray} \tag{ 22 }$

where ${\tau }_{\mathrm{IR}}$ is the infrared optical depth at the midplane, and

$\begin{eqnarray}{T}_{{\rm{s}},\mathrm{vis}} & = & {\left(\displaystyle \frac{3{{GM}}_{* }\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{R}^{3}}\right)}^{1/4}\\ & \approx & 107{\left(\displaystyle \frac{R}{1\mathrm{au}}\right)}^{-3/4}{\left(\displaystyle \frac{{M}_{* }}{2.5{M}_{\odot }}\right)}^{1/4}{\left(\displaystyle \frac{\dot{M}}{{10}^{-8}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{1/4}\,{\rm{K}}\end{eqnarray} \tag{ 23 }$

is the surface temperature determined by the mass accretion rate $\dot{M}=3\pi {{\rm{\Sigma }}}_{{\rm{g}}}\nu$ (e.g., Shakura & Sunyaev 1973), where ν is the turbulent viscosity.

To see when viscous heating can be neglected, we compare Equation (22) with (15), which is the minimum estimate for the midplane temperature of a passively irradiated disk (see Figure 7). Combining Equations (15) and (22), together with ${\tau }_{\mathrm{IR}}=(1/2){\kappa }_{{\rm{d}}}{{\rm{\Sigma }}}_{{\rm{d}}}=(1/2){f}_{{\rm{d}}2{\rm{g}}}{\kappa }_{{\rm{d}}}{{\rm{\Sigma }}}_{{\rm{g}}}$ (where ${{\rm{\Sigma }}}_{{\rm{g}}}$ and ${{\rm{\Sigma }}}_{{\rm{d}}}$ are the surface densities of gas and dust), we find that ${T}_{\mathrm{mid},{\rm{D}}}\gt {T}_{\mathrm{mid},\mathrm{vis}}$ if

$\begin{eqnarray}&&0.8{\left(\displaystyle \frac{{L}_{* }}{56{L}_{\odot }}\right)}^{-2/7}{\left(\displaystyle \frac{{M}_{* }}{2.5{M}_{\odot }}\right)}^{31/70}{\left(\displaystyle \frac{{\kappa }_{{\rm{d}}}}{700{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{1/5}{\left(\displaystyle \frac{{f}_{{\rm{d}}2{\rm{g}}}}{{10}^{-2}}\right)}^{1/5}\\ &&\quad \times {\left(\displaystyle \frac{\dot{M}}{{10}^{-8}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{2/5}{\left(\displaystyle \frac{\alpha }{{10}^{-3}}\right)}^{-1/5}{\left(\displaystyle \frac{R}{1\mathrm{au}}\right)}^{-33/70}\lt 1.\end{eqnarray} \tag{ 24 }$

Here, we have used the α-prescription for the turbulent viscosity (Shakura & Sunyaev 1973), i.e., $\nu =\alpha {c}_{{\rm{s}}}{h}_{{\rm{g}}}$ . We have also assumed ${\tau }_{\mathrm{IR}}\gg 1$ , ${z}_{* ,{\rm{D}}}/{h}_{{\rm{g}}}=3.6$ , and ${T}_{{\rm{s}},{\rm{D}}}^{4}=2{T}_{2}^{4}$ for simplicity. Equation (24) suggests that in a disk around a Herbig Ae star of ${L}_{* }=56{L}_{\odot }$ and ${M}_{* }=2.5{M}_{\odot }$ , and $\dot{M}={10}^{-8}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ , stellar irradiation dominates over viscous heating if ${f}_{{\rm{d}}2{\rm{g}}}\lesssim {10}^{-2}$ , in agreement with the result of the hydrodynamical simulation by Flock et al. (2016; see their Figure 8). Even in an inner part $(\sim 0.1\,\mathrm{au})$ of a disk around a T-Tauri star of ${L}_{* }={L}_{\odot }$ and ${M}_{* }={M}_{\odot }$ , stellar irradiation dominates if ${f}_{{\rm{d}}2{\rm{g}}}\lesssim {10}^{-3}$ and $\dot{M}\lesssim 6\times {10}^{-10}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ . Therefore, our temperature model is applicable to weakly accreting T-Tauri disks where tiny, opacity-dominating dust is depleted (through, e.g., dust coagulation).

3.3. Migration Trap

A steep temperature drop as seen in the borders of regions B, C, and D is known to act as a planet trap where the strong corotation torque halts the inward migration of planets due to the negative Lindblad torque from the gas disk (e.g., Paardekooper et al. 2010). In addition, if the boundary between regions C and D sets the dead zone inner edge as discussed in Section 3.1, the positive surface density gradient arising from the change in turbulent viscosity also prevents planets from inward migration (Masset et al. 2006).

It is meaningful to compare the locations of the borders to the orbital distribution of observed planets in order to reveal whether the borders really act as the migration trap. Our solutions show that the positions of the boundaries between regions B, C, and D are roughly proportional to ${L}_{* }^{1/2}$ , though the position of the boundary between regions C and D has the additional factor $\sqrt{1+{\rm{\Gamma }}}$ , which slightly depends on the stellar parameters. If we assume ${L}_{* }\propto {M}_{* }^{3/2}$ (Mulders et al. 2015), these boundaries are scaled as ${M}_{* }^{3/4}$ . Mulders et al. (2015) investigated Kepler Objects of Interest (KOIs) and showed that the distance from the star where the planet occurrence rate drops scales with semimajor axis as ${M}_{* }^{1/3}$ , which is different from our estimate. However, the spectral range is not enough for a detailed comparison because the majority of the host stars for KOIs are F, G, and K stars. In order to discuss how the boundaries really act as a migration trap, observations of planets around A and B stars are needed. Mulders et al. (2015) also showed that many KOIs orbit at ∼0.1 au ,which is similar to the value of ${R}_{\mathrm{DIB}}$ estimated from our model.

3.4. Instabilities at the Inner Region of Protoplanetary Disks

Steep temperature gradients at the boundaries between regions B, C, and D would induce various kind of instabilities. One possible instability is the subcritical baroclinic instability, which is the instability powered by radial buoyancy due to the unstable entropy gradient (Klahr & Bodenheimer 2003; Lyra 2014). If we assume the power-law density and temperature profile as $\rho \propto {R}^{-{\beta }_{\rho }}$ and $T\propto {R}^{-{\beta }_{T}}$ , respectively, the condition for the subcritical baroclinic instability is given by (Klahr & Bodenheimer 2003)

$\begin{eqnarray}&&{\beta }_{T}-({\gamma }_{\mathrm{ad}}-1){\beta }_{\rho }\gt 0,\end{eqnarray} \tag{ 25 }$

where ${\gamma }_{\mathrm{ad}}$ is the adiabatic index. This condition would be satisfied at the boundaries between regions B, C, and D, where ${\beta }_{T}$ is large. The boundary between regions C and D might also have a positive surface density gradient ( ${\beta }_{\rho }\lt 0$ ) if the dead zone edge inner edge lies there. The positive density gradient would further enhance the instability on the boundary. However, the subcritical baroclinic instability has not been observed in Flock et al. (2016). This is because the spatial resolution of Flock et al. (2016) is not enough to resolve the instability. Lyra (2014) showed that for resolutions $\gt 32$ cells/ ${h}_{{\rm{g}}}$ the subcritical instability converges, while Flock et al. (2016) used around 15 cells/ ${h}_{{\rm{g}}}$ . Also, Flock et al. (2016) calculated for around 20 local orbits at $2\,\mathrm{au}$ , which is not enough for the subcritical instability to operate. The growth timescale for the instability requires many hundreds of local orbits (Lyra 2014).

The Rossby wave instability could also be induced at the boundary between regions C and D due to the steep surface density gradient produced by the dead zone edge (Lyra & Mac Low 2012). In fact, the magnetohydrodynamical simulation by Flock et al. (2017) shows that the boundary between regions C and D produces a vortex, indicating that the Rossby wave instability should operate there.

If such instabilities operate on these boundaries, vortices caused by the instabilities would accumulate dust (e.g., Barge & Sommeria 1995), which in turn could lead to rocky planetesimal formation via the streaming instability (e.g., Youdin & Goodman 2005). We will explore this possibility in future work.

4. Summary

We have analytically derived the temperature and dust-to-gas mass ratio profile for the inner region of protoplanetary disks based on the results from the recent hydrodynamical simulations conducted by Flock et al. (2016). The temperature profile for the inner region of protoplanetary disks can be divided into four regions. The innermost region is dust-free and optically thin, with the temperature determined by the gas opacity (Equation (3)). As the temperature goes down and approaches the dust evaporation temperature, silicate dust starts to condense, producing an optically thin dust halo with a nearly constant temperature regulated by partial dust condensation. We have derived the dust-to-gas mass ratio profile in the dust halo using the fact that partial dust condensation regulates the temperature to the dust evaporation temperature (Equation (4)). Beyond the dust halo, there is an optically thick condensation front where all the available silicate gas condenses out. The curvature of the condensation surface is simply determined by the condition that the surface temperature must be nearly equal to the characteristic temperature ∼1200 K (Equation (9)). The temperature profile for the outermost region is essentially the same as the classical optically thick temperature profile (e.g., Kusaka et al. 1970). We have derived the midplane temperature in the outer two regions using the two-layer approximation, with additional heating by the condensation front for the outermost region. As a result, the overall temperature profile follows a step-like profile with steep temperature gradients at the borders between the outer three regions. The borders might act as planet traps where the inward migration of planets due to gravitational interaction with the gas disk stops. The temperature at the border between the two outermost regions coincides with the temperature needed to activate magnetorotational instability, suggesting that the inner edge of the dead zone must lie at this border. The radius of the dead zone inner edge predicted from our temperature profile is ∼2–3 times larger than that expected from the classical optically thick temperature.

We thank Neal Turner and Shigeru Ida for insightful discussions. We also thank the anonymous referee for useful comments on our manuscript. S.O. is supported by Grants-in-Aid for Scientific Research (No. 15H02065, 16H04081, 16K17661).

Appendix: Appendix Effect of the Evaporation Temperature

Dust evaporation temperature ${T}_{\mathrm{ev}}$ is modeled as a function of the gas density ρ as (Isella & Natta 2005)

$\begin{eqnarray}&&{T}_{\mathrm{ev}}=2000{\left(\displaystyle \frac{\rho }{1{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{0.0195}\,{\rm{K}}.\end{eqnarray} \tag{ 26 }$

Using the relation $\rho ={{\rm{\Sigma }}}_{{\rm{g}}}/\sqrt{\pi }{h}_{{\rm{g}}}$ and substituting Equation (26) as ${T}_{{\rm{B}}}$ , Equations (5) and (6) can be rewritten as

$\begin{eqnarray}{R}_{\mathrm{AB}} & = & \left[0.079\left(\displaystyle \frac{{R}_{* }}{0.01\,\mathrm{au}}\right){\left(\displaystyle \frac{{T}_{* }}{5000{\rm{K}}}\right)}^{2}\right.\\ & & \times {\left.{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{100{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{-0.0386}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-0.0193}\right]}^{\gamma }\,\mathrm{au}\end{eqnarray} \tag{ 27 }$

and

$\begin{eqnarray}&&{R}_{\mathrm{BC}}={\left(\displaystyle \frac{{\kappa }_{{\rm{d}}}({T}_{* })}{{\kappa }_{{\rm{d}}}(T)}\right)}^{\gamma /2}{R}_{\mathrm{AB}},\end{eqnarray} \tag{ 28 }$

respectively. Here, ${{\rm{\Sigma }}}_{0}$ is the surface density at 1 au. The index γ is related to the index of the radial surface density profile, ${{\rm{\Sigma }}}_{{\rm{g}}}\propto {r}^{-\beta }$ , as $\gamma =2.0195/(1.9025-0.078\beta )$ . For $\beta =1$ , $\gamma \approx 1.107$ . Equation (27) is similar to the expression derived by (Kama et al. 2009, their Equation (A.4)), although the β dependence is different. If the disk is massive ( ${{\rm{\Sigma }}}_{0}\gg 1000\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ ), the boundary between regions B and C (i.e., the radial position for ${\tau }_{* }=1$ ) moves closer to the star compared to ${R}_{\mathrm{BC}}$ written by Equations (6) and (28) because of the large optical depth. In this case, in order to determine the radial position of ${\tau }_{* }=1$ , we have to calculate the optical depth using Equation (4).