# Radiation Hydrodynamics Models of the Inner Rim in Protoplanetary Disks

M. Flock, S. Fromang, N. J. Turner, M. Benisty

► To cite this version:

M. Flock, S. Fromang, N. J. Turner, M. Benisty. Radiation Hydrodynamics Models of the Inner Rim in Protoplanetary Disks. The Astrophysical Journal, 2016, 827, 10.3847/0004-637X/827/2/144. insu-03691528

HAL Id: insu-03691528

https://insu.hal.science/insu-03691528v1

Submitted on 21 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RADIATION HYDRODYNAMICS MODELS OF THE INNER RIM IN PROTOPLANETARY DISKS

M. FLOCK¹,², S. FROMANG², N. J. TURNER¹, AND M. BENISTY³

¹Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA; mflock@caltech.edu

²Laboratoire AIM, CEA/DSM-CNRS-Université Paris 7, Irfu/Service d'Astrophysique, CEA-Saclay, F-91191 Gif-sur-Yvette, France

³Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

Received 2016 February 25; revised 2016 April 15; accepted 2016 April 15; published 2016 August 18

ABSTRACT

Many stars host planets orbiting within a few astronomical units (AU). The occurrence rate and distributions of masses and orbits vary greatly with the host star's mass. These close planets' origins are a mystery that motivates investigating protoplanetary disks' central regions. A key factor governing the conditions near the star is the silicate sublimation front, which largely determines where the starlight is absorbed, and which is often called the inner rim. We present the first radiation hydrodynamical modeling of the sublimation front in the disks around the young intermediate-mass stars called Herbig Ae stars. The models are axisymmetric and include starlight heating; silicate grains sublimating and condensing to equilibrium at the local, time-dependent temperature and density; and accretion stresses parameterizing the results of MHD magnetorotational turbulence models. The results compare well with radiation hydrostatic solutions and prove to be dynamically stable. Passing the model disks into Monte Carlo radiative transfer calculations, we show that the models satisfy observational constraints on the inner rim's location. A small optically thin halo of hot dust naturally arises between the inner rim and the star. The inner rim has a substantial radial extent, corresponding to several disk scale heights. While the front's overall position varies with the stellar luminosity, its radial extent depends on the mass accretion rate. A pressure maximum develops near the location of thermal ionization at temperatures of about 1000 K. The pressure maximum is capable of halting solid pebbles' radial drift and concentrating them in a zone where temperatures are sufficiently high for annealing to form crystalline silicates.

**Key words:** accretion, accretion disks – magnetohydrodynamics (MHD) – radiative transfer – techniques: interferometric

1. INTRODUCTION

Thousands of the planets discovered with *Kepler* and ground-based telescopes orbit within an astronomical unit of their low-mass host stars (Johnson et al. 2013; Benz et al. 2014, p. 691; Lissauer et al. 2014). In contrast, intermediate-mass stars around 2 solar masses more often host super-Jovian planets orbiting outside 1 AU (Bowler et al. 2010), while stars of more than 2.7 solar masses have few or no super-Jovians at these distances (Reffert et al. 2015). Our understanding of these diverse planetary systems' origins relies on our knowledge of the central regions of the disks around all young stars.

One key location in the disks is the silicate sublimation front, the boundary between transparent rock vapor and an opaque dust cloud, and thus between hot gas lit directly by the star and warm material in the shadow (Dullemond & Monnier 2010).

Another key location, lying near the sublimation front but distinct from it, is the turbulent front. This surface separates material ionized well enough to undergo magnetorotational instability (MRI), leading to turbulence, from neutral material that is laminar or subject to weak turbulence driven by hydrodynamical instabilities (Balbus & Hawley 1998; Turner et al. 2014). Disk annuli experiencing magnetorotational turbulence have lower surface densities than their less turbulent neighbors if the overall flow is in steady state, since the MRI-active regions' stronger accretion stresses drive the gas through faster. Thus, the turbulent front leads to a nearby local maximum in the radial profile of the midplane gas pressure (Varnière and Tagger 2006; Dzyurkevich et al. 2010).

The pressure maximum collects solid particles, because gas outside the peak orbits slower than Keplerian, giving the particles a headwind so they drift inward, while gas inside the

peak orbits faster than Keplerian, yielding a tailwind that raises particles' orbits (Weidenschilling 1977; Haghighipour & Boss 2003; Lyra et al. 2008, 2009). In particular, the turbulent front can affect the distribution of planet-forming solids in intermediate-mass Herbig stars' disks (Kretke et al. 2009). The turbulent front's surface density jump also provides favorable conditions for the growth of cyclonic vortices, which can further concentrate solid particles along the azimuthal direction (Barge & Sommeria 1995; Lyra & Mac Low 2012; Faure et al. 2014). Concentrating the solid material makes grain growth efficient (Testi et al. 2014, p. 339), and it could lead to in situ planet formation (Chatterjee & Tan 2014), especially at the inner disk (Boley et al. 2014). Furthermore, the pressure maximum can halt young planets' orbital migration (Masset et al. 2006; Matsumura et al. 2009; Lyra et al. 2010; Kretke & Lin 2012; Bitsch et al. 2014; Hu et al. 2015). The turbulent front thus plays several roles in planetary systems' development.

The sublimation and turbulent fronts are worth considering together when modeling the planet-forming environments near young stars, because they are mutually coupled. The sublimation front affects the turbulent front by governing the starlight absorption and heating. The turbulent front in turn affects the sublimation front through its control over the radial distribution of material.

Interferometric observations of Herbig Ae/Be stars at near-infrared (NIR) wavelengths can resolve the sublimation and turbulent fronts' locations (Dullemond & Monnier 2010; Kraus 2015). However, interpreting the measured surface brightnesses in terms of the disk's density and temperature structure remains a challenge. The ingredients that must be
considered include the transfer of the starlight and the infrared photons re-radiated by the disk, the dust particles' sublimation and condensation, and the gas vertical hydrostatic equilibrium (Kama et al. 2009). Early radiation hydrostatic models had difficulty fitting Herbig stars' spectral energy distributions (SEDs), especially underpredicting the flux at NIR wavelengths (Hillenbrand et al. 1992; Chiang et al. 2001; Meeus et al. 2001; Millan-Gabet et al. 2001). Attempts to solve this problem began with a vertical wall of dust (Dullemond et al. 2001; Natta et al. 2001). The vapor between wall and star was optically thin, letting the starlight fall directly on the wall. This made the wall extra hot and thus extra tall under hydrostatic equilibrium. The puffed-up wall intercepted extra starlight, yielding a higher NIR flux. The next generation of models included the sublimation temperature's density dependence, leading to a rounded shape for the sublimation front: the wall of dust sloped away from the star above and below the midplane, improving the match to disks inclined at a range of angles with respect to our line of sight (Isella & Natta 2005; Isella et al. 2006).

Most works have neglected the absorption of starlight in the hot gas interior to the sublimation front. However, Muzerolle et al. (2004) suggested that gas between the sublimation front and the star can be optically thick enough to push the sublimation front inward significantly. Gas interior to the sublimation front in the disks around several Herbig Ae/Be stars was recently detected by CO rovibrational emission (Ilee et al. 2014) or Brγ emission (Mendigutía et al. 2015). An issue not yet addressed by any model is whether the sublimation front is dynamically stable. Does the sharp temperature jump destabilize displacements at the front? In this work we investigate the structure and stability of the sublimation and turbulent fronts, using radiation hydrostatic and, for the first time, radiation hydrodynamical models. The models are axisymmetric and treat the transfer of both starlight and infrared radiation, with the opacity linked to the grains' sublimation and condensation. The accretion stresses and heating come from a prescription for MRI turbulence with a switch at the onset of thermal ionization. All the models have surface density profiles corresponding to steady inflow past the two fronts toward the star. We first consider radiation hydrostatic models, then radiation hydrodynamical models. Finally, we post-process the results through Monte Carlo radiative transfer calculations to construct synthetic observations of the disk model, which we compare with constraints from observations.

The structure of the paper is as follows. In Section 2 we present the radiation hydrostatic method, the dust and gas opacity, and the dust sublimation module. In Section 3 we present a first-principle 2D radiation hydrostatic solution to define the inner rim and to show the temperature and density structure, followed by a steady-state model in Section 4. In Section 5 we investigate the effects of the model parameter, and finally we present the radiation hydrodynamical models in Section 6. In Section 7, we compare our results with observational constraints and calculate the SEDs of our models, followed by a discussion (Section 8) and our conclusions (Section 9).

## 2. METHOD

The aim of the method can be stated as the following: given a star with known properties (mass, radius, luminosity) and a steady-state mass accretion rate $\dot{M}$ onto that star, what are the

spatial distributions of matter (both gas and dust) and temperature?

Here, the three important timescales are the radiative timescale, $t_{\text{rad}}$, over which the temperature reaches equilibrium (given a dust and gas density field); a dynamical timescale associated with sound wave propagation, $t_{\text{dyn}}$, over which the vertical disk equilibrium is set; and a long timescale, $t_{\text{visc}}$, associated with angular momentum transport in the disk (and often referred to as the viscous timescale), over which the disk surface density $\Sigma$ evolves. The method we use in the present paper relies on the large decoupling between the different timescales: $t_{\text{rad}}^{\tau=1} \ll t_{\text{dyn}} \ll t_{\text{rad}}^{\tau \gg 1} \ll t_{\text{visc}}$, given the optical thickness $\tau$ of the thermal emission in the disk.

Because of this ordering, the local volume densities remain near vertical hydrostatic balance while the total surface density profile evolves. As we are interested in the steady-state surface density, we determine $\Sigma$ by assuming a uniform mass accretion rate $\dot{M}$. From the surface density, we calculate the dust and gas densities jointly with the temperature, using an iterative method.

In Section 2.1 we summarize the basic iterative procedure described by Flock et al. (2013). In the present paper, we modify this procedure in two key ways. First, because we seek to model the sublimation front, we let the dust-to-gas ratio $f_{\text{D2G}}$ vary with position, greatly affecting the dust opacities. Second, the sublimation front lies near the thermal ionization front, across which magnetic stresses can vary abruptly. Over the viscous timescale, such a jump in the stress ought to determine the surface density profile. We describe the procedure for computing $f_{\text{D2G}}$ in Section 2.2 and the procedure for estimating the steady-state surface density profile $\Sigma$ in Section 2.3. Finally, we outline how we compute the dust and gas opacities in Section 2.4.

### 2.1. Hydrostatic Disk Structure

The iterative method solves the coupled equations that describe vertical hydrostatic balance at a given spherical radius $r$ in the disk. In a spherical coordinate system $(r, \theta, \phi)$, these equations are

$$ \frac{\partial P}{\partial r} = -\rho \frac{\partial \Phi}{\partial r} + \frac{\rho v_{\phi}^2}{r} \qquad (1) $$

$$ \frac{1}{r} \frac{\partial P}{\partial \theta} = \frac{1}{\tan \theta} \frac{\rho v_{\phi}^2}{r}, \qquad (2) $$

where $\rho$ is the gas density, $v_{\phi}$ is the gas azimuthal velocity, the gravitational potential $\Phi = GM_{*}/r$ with the gravitational constant $G$, stellar mass $M_{*}$, and $P$ is the thermal pressure that relates to the temperature $T$ through the ideal gas equation of state:

$$ P = \frac{\rho k_{\text{B}} T}{\mu_{\text{g}} u}, \qquad (3) $$

with the mean molecular weight $\mu_{\text{g}}$, the Boltzmann constant $k_{\text{B}}$, and the atomic mass unit $u$. The temperature in the disk is set by a balance between stellar irradiation and radiative cooling. For a given density field, the radiation equilibrium is obtained as the steady-state solution to the following coupled
set of equations:

$$ \frac{1}{\Gamma - 1} \partial_t P = -\sigma_p c (a_R T^4 - E_R) - \nabla \cdot F_*, \\ \partial_t E_R - \nabla \frac{c\lambda}{\sigma_R} \nabla E_R = +\sigma_p c (a_R T^4 - E_R), \qquad (4) $$

with adiabatic index $\Gamma$, radiation energy $E_R$, irradiation flux $F_*$, Rosseland and Planck opacity $\sigma_R$ and $\sigma_p$, radiation constant $a_R = 4\sigma_b/c$ with Stefan–Boltzmann constant $\sigma_b$, and speed of light $c$. The flux limiter

$$ \lambda = \frac{2 + R_{\lambda}}{6 + 3R_{\lambda} + R_{\lambda}^2} \qquad (5) $$

is taken from Levermore & Pomraning (1981, Equation (28) therein) with

$$ R_{\lambda} = \frac{|\nabla E_{R}|}{\sigma_{R} E_{R}}. \qquad (6) $$

The gas is a mixture of molecular hydrogen and helium with solar abundance (Decampli et al. 1978) so that $\mu_g = 2.35$ and $\Gamma = 1.42$. In Appendix A, we provide a list of all the constants we use.

In this work we consider the frequency-integrated irradiation flux $F_*$ at radius $r$:

$$ F_{*}(r) = \left(\frac{R_{*}}{r}\right)^{2} \sigma_{b} T_{*}^{4} e^{-\tau_{*}}, \qquad (7) $$

with $T_*$ and $R_*$ being the surface temperature and the radius of the star, respectively. The radial optical depth of the irradiation flux is defined at each position $\theta$ by

$$ \tau_{*}(r) = \int_{R_{*}}^{r} \sigma_{*} dr = \tau_{0} + \int_{r_{0}}^{r} \sigma_{*} dr, \qquad (8) $$

where $r_0$ denotes the computational box inner radius and $\sigma_*$ is the opacity at the stellar temperature (see below). $\tau_0$ is the inner optical depth located between the surface of the star and $r_0$ and determines how much of the incoming irradiation is blocked before entering the computational domain. We assume a pure gas disk located between three stellar radii and $r = r_0$, which gives $\tau_0 = \kappa_{\text{gas}} \rho_{r_0} (r_0 - 3R_*)$. $\sigma_*$ is determined by $\sigma_* = \rho_{\text{dust}} \kappa_P (T_*) + \rho_{\text{gas}} \kappa_{\text{gas}}$. The dust Planck opacity $\kappa_P$ and the gas opacity are defined in Section 2.4. The inner gas disk edge is assumed to be located at 3 stellar radii, which is the position of the magnetospheric truncation radius for a Herbig-type star (Muzerolle et al. 2004). We note that this ray-tracing approach neglects scattered starlight.

The iterative method is summarized below: first, we set the surface density profile $\Sigma(R)$ at the cylindrical radius $R$, and a temperature field $T(r, \theta)$, which is calculated using the optically thin solution. We then calculate $\rho(r, \theta)$ and $v_\phi(r, \theta)$ by solving Equations (1) and (2). We next get the new temperature profile from the new radiation equilibrium solution of Equation (4). We iterate the last two steps until we reach convergence. The reader is referred to Section 3.1 of Flock et al. (2013) for more details.

## 2.2. Dust Sublimation

The inner disk’s structure depends critically on when and where the dust sublimates. Sublimation also complicates obtaining converged solutions with iterative methods (Kama et al. 2009; Vinković 2012) for at least two reasons. First, some material is heated by radiation from outside as well as from inside its orbit. In particular, grains starward from the front are lit on their night sides by the infrared radiation from their more distant neighbors. This “backwarming” leaves them hotter than they would be in optically thin surroundings, and this effect is self-consistently included in our method. Second, the front is geometrically very thin, since a small column of grains suffices to shield the material beyond from the starlight, making the transition from vapor to condensed quite sharp. High spatial resolution is needed to resolve the layer where the stellar flux is deposited. In our method, we smooth this transition, which allows us to overcome this resolution constraint as the models should be suitable for future fully 3D radiation MHD simulations.

We follow Pollack et al. (1994) and the fitting model of Isella & Natta (2005)<fnref n="4" /> that applies to situations for which the most refractory grains are silicates. Then the dust sublimation temperature $T_{\text{ev}}$ is set by

$$ T_{\text{ev}} = 2000 \text{ K} \left( \frac{\rho}{1 \text{ g cm}^{-3}} \right)^{0.0195} \qquad (9) $$

$T_{\text{ev}}$ is then used to calculate the dust-to-gas ratio $f_{\text{D2G}}$, i.e., the ratio between the dust density and the gas density. We use

$$ f_{\text{D2G}} = \left\{ \frac{f_{\Delta\tau}}{2} \left\{ 1 - \tanh \left[ \left( \frac{T - T_{\text{ev}}}{100 \text{ K}} \right)^3 \right] \right\} \left\{ \frac{1 - \tanh(1 - \tau_{*})}{2} \right\} \right. \\ \text{if } T > T_{\text{ev}} \quad f_0 \left\{ \frac{1 - \tanh(20 - \tau_{*})}{2} \right\} + f_{\Delta\tau} \\ \text{if } T < T_{\text{ev}}, \qquad (10) $$

with the reference dust-to-gas mass ratio $f_0 = 0.01$ and $f_{\Delta\tau} = 0.3/(\rho_{\text{gas}} \kappa_P \Delta r)$ setting the dust amount to account for an optical depth of $\Delta\tau_* = 0.3$. Such a value ensures the resolution of the absorption of the irradiation at the rim. For $T > T_{\text{ev}}$, Equation (10) is similar to the formula of Kama et al. (2009) and controls the amount of dust for temperatures above the sublimation temperature. Note that, in this regime, we also impose a minimum value of the dust-to-gas mass ratio: $f_{\text{D2G}}^{\min} = 10^{-10}$. The upper limit, in this regime, is given by the value $f_{\Delta\tau}$, which is reached close to $T_{\text{ev}}$ and $\tau_* = 1$. For $T < T_{\text{ev}}$, Equation (10) limits the dust amount until the irradiation is absorbed, which is reached close to $\tau_* = 20$. Finally, the dust-to-gas mass ratio computed by $f_{\text{D2G}}$ is then used to define the total opacity in each grid cell for the irradiation and for the thermal emission (see Section 2.4). We note that the opacity gradient across the rim could be even more gradual than we model here if the dust consists of components with differing sublimation thresholds (McClure et al. 2013). In addition, the species likely to occur in protostellar disks cover a significant range in sublimation temperature (Pollack et al. 1994).

<footnote n="4">⁴ We note that a slightly different fit was presented by Vinković (2014) (see Equation (17) therein). We have tested both expressions and found no significant difference on the rim structure.</footnote>
In Appendix B, we show the robustness of this function by comparing different dust sublimation functions for $T > T_{ev}$. We demonstrate the importance of the function for $T < T_{ev}$ in Appendix C. A resolution study is performed in Appendix F.

## 2.3. Surface Density Radial Profile

The disk surface density radial profile is governed by the transport of angular momentum and evolves over long timescales of several thousands of orbits. Here our strategy is to use simple scaling laws that come from our understanding of the angular momentum transport to obtain a crude estimate of the equilibrium $\Sigma(R)$ and use it along with the iterative method described above to compute the disk structure.

The angular momentum transport in protoplanetary disks' inner regions is most likely dominated by MHD turbulence mediated by the MRI. As a result, the flow is very complex and time dependent. A simple, yet efficient, way to describe this complex flow is the $\alpha$-prescription, which mimics the effect of the turbulence with a viscosity $\nu_t$ (Shakura & Sunyaev 1973). It is customary to write this “turbulent” viscosity as

$$ \nu_t = \frac{\alpha c_s^2}{\Omega}, \qquad (11) $$

with local sound speed $c_s$ and disk rotation frequency $\Omega = \sqrt{GM_*/R^3}$. Using this large-scale model in steady state and far away from the central star, $\dot{M}$ and $\nu_t$ are related by

$$ \Sigma(R) = \frac{\dot{M}}{3\pi\nu_t(R)}. \qquad (12) $$

For a given value of $\dot{M}$, $\Sigma$ is thus smaller in the turbulent regions of the disk (where $\alpha$ is large) than in the laminar parts (where $\alpha$ is small). In this paper, we use Equation (12) to estimate $\Sigma(R)$. The viscosity is determined with Equation (11), in which we specify $\alpha$ using

$$ \alpha = (\alpha_{\text{in}} - \alpha_{\text{out}}) \left[ \frac{1 - \tanh\left(\frac{T_{\text{MRI}} - T}{\Delta T}\right)}{2} \right] + \alpha_{\text{out}}, \qquad (13) $$

where $T$ here stands for the midplane temperature of the disk. This formula ensures that $\alpha$ varies smoothly from $\alpha_{\text{in}}$ at those locations where the MRI is active ($T > T_{\text{MRI}}$) to $\alpha_{\text{out}}$ for $T < T_{\text{MRI}}$. The threshold temperature $T_{\text{MRI}}$ for magnetorotational turbulence is about 1000 K, while the exact value depends on the dust-to-gas ratio, the grain sizes, and the gas density (Umebayashi & Nakano 1988; Desch & Turner 2015). The likely ranges in these parameters allow threshold temperatures between about 800 and 1200 K, and we study the effect of these values of $T_{\text{MRI}}$ on the results in Section 5.4.

In the rest of this paper, we use $T_{\text{MRI}} = 1000$ K, $\alpha_{\text{in}} = 1.9 \times 10^{-2}$, $\alpha_{\text{out}} = 10^{-3}$, and $\Delta T = 25$ K. By varying $\Delta T$ from 10 to 50 K, we have checked that the influence of its exact value on the results is small and does not modify our conclusions. A value of $\alpha_{\text{in}}$ slightly above 0.01 is justified by MRI simulations in well-ionized media (Fromang & Nelson 2006; Davis et al. 2010; Simon et al. 2011; Flock et al. 2012). However, the value of $\alpha_{\text{out}}$ is poorly known and depends on the dominant nonideal MHD term (Lesur et al. 2014; Turner et al. 2014; Simon

et al. 2015) and the strength of possible hydrodynamic drivers of turbulence (Nelson et al. 2013; Klahr & Hubbard 2014; Lyra 2014). We study the sensitivity of the results to the precise value of $\alpha_{\text{out}}$ in Section 5.3.

## 2.4. Opacities

In this paper, we consider gas and dust opacity. We assume $\sigma_R = \sigma_P = \rho_{\text{dust}}\kappa_P(T) + \rho_{\text{gas}}\kappa_{\text{gas}}$, and we simplify the problem by reducing the frequency-dependent opacities to three different frequency-averaged values: the gas opacity $\kappa_{\text{gas}}$, the Planck dust opacity at the rim $\kappa_P(T_{\text{rim}})$, and the Planck dust opacity at the stellar temperature $\kappa_P(T_*)$. Such a simplified model compares very well with a more complex model including frequency-dependent irradiation and temperature-dependent dust opacity (see Appendix D for the full comparison). In the following we detail our choices for the dust and gas opacities.

### 2.4.1. Dust Opacity

We generate the dust opacity table for different wavelengths using the MieX code by Wolf & Voshchinnikov (2004). For details on the dust material, the opacity calculations, and comparison to other dust opacity tables, we refer to Appendix A. As mentioned above, two wavelength bands and thus two opacity values are important at the rim. The first is the dust opacity at the irradiation temperature. For $T_* = 10,000$ K we set $\kappa_P(T_*) = 2100 \text{ cm}^2 \text{ g}^{-1}$ (exact value $\kappa_P(10,000 \text{ K}) = 2100.3 \text{ cm}^2 \text{ g}^{-1}$), which is the mean opacity of the dust for the stellar spectrum. We note that we will consider different stellar types, and so the value of $\kappa_P(T_*)$ will be adapted according to the surface temperature of the star. The second frequency-averaged opacity is the dust opacity at the dust sublimation temperature. Typical sublimation temperatures are between 1300 K ($\kappa_P(1300 \text{ K}) = 690.1 \text{ cm}^2 \text{ g}^{-1}$) and 1400 K ($\kappa_P(1400 \text{ K}) = 717.9 \text{ cm}^2 \text{ g}^{-1}$). We fix the second dust opacity, which represents the cooling efficiency at the rim, to a value of $\kappa_P(T_{\text{rim}}) = 700 \text{ cm}^2 \text{ g}^{-1}$.

In this context, the ratio between emission and absorption efficiency $\epsilon$ is important as it determines the dust temperature in an optically thin, stellar irradiated environment. The small dust particles have a high opacity at short wavelengths compared to longer wavelengths. They are more efficient in absorbing the shorter wavelength of the stellar radiation than in emitting at longer thermal wavelengths, so they appear hotter than a perfect blackbody radiator. For our model, the ratio of the emission-to-absorption efficiency of the dust is then given by $\epsilon = \kappa_P(T_{\text{rim}})/\kappa_P(T_*) = 1/3$. Such a value is typical for a mixture of dust particles with different sizes (Dullemond & Monnier 2010). Small particles of single size 0.1 µm have a value of $\epsilon = 0.08$ (Isella & Natta 2005; Kama et al. 2009), while including larger particles increases the value of $\epsilon$.

### 2.4.2. Gas Opacity

The gas opacity is more difficult to compute as it is dominated by the various line opacities (Helling et al. 2000; Dullemond & Monnier 2010). A fine frequency sampling is required to recover correct values for the mean Planck opacities (Malygin et al. 2014). In addition, such mean opacities can
<table><caption>**Table 1**

General Setup Parameter for the Radiation Hydrostatic Disk Model S100</caption><tr><td>Surface density</td><td>100 g cm<sup>-2</sup>, uniform</td></tr><tr><td>N<sub>r</sub> × N<sub>θ</sub></td><td>1280 × 128</td></tr><tr><td>Cell aspect ratio</td><td>rΔθ/Δr ~ 1.2</td></tr><tr><td>r<sub>in</sub> - r<sub>out</sub> : Z/R</td><td>0.2–4 AU : ~±0.18</td></tr><tr><td>Stellar parameter</td><td>T<sub>*</sub> = 10000 K, R<sub>*</sub> = 2.5 R<sub>☉</sub><br>M<sub>*</sub> = 2.5 M<sub>☉</sub></td></tr><tr><td>Opacity</td><td>κ<sub>P</sub>(T<sub>*</sub>) = 2100 cm<sup>2</sup> g<sup>-1</sup><br>κ<sub>P</sub>(T<sub>rim</sub>) = 700 cm<sup>2</sup> g<sup>-1</sup><br>κ<sub>gas</sub> = 10<sup>-4</sup> cm<sup>2</sup> g<sup>-1</sup></td></tr><tr><td>Dust-to-gas mass ratio</td><td>f<sub>0</sub> = 0.01</td></tr></table>

become very high (see Figure 2 in Malygin et al. 2014). This is because the gas opacity at smaller wavelengths ($\lambda < 0.45$ µm) is high due to H$_2$ and CO transitions for wavelengths shorter than ~0.45 µm (Muzerolle et al. 2004). However, as Muzerolle et al. (2004) pointed out, there is a lack of absorption in a wavelength range for $\lambda > 0.45$ µm for the stellar irradiation. Dependent on the molecular abundance and composition, the gas opacity can vary between $10^{-6}$ and $1$ cm$^2$ g$^{-1}$ (Dullemond & Monnier 2010; Malygin et al. 2014) for NIR wavelengths. In this work, we fix the gas opacity to $\kappa_{gas} = 10^{-4}$ cm$^2$ g$^{-1}$. The value is chosen so that the radial optical depth $\tau_{*}$ and vertical optical depth $\tau_{Z}^{\text{NIR}}$ remain small, with

$$ \tau_{Z}^{\text{NIR}} = \int_{-\infty}^{\infty} \sigma_{p} dz. \qquad (14) $$

This is very important as otherwise the inner gas disk would block the irradiation and the rim radius would move too close to the star (see also Section 5.1), inconsistent with observations. We note also that we use the same gas opacity value for the irradiation and for the thermal emission, which results in the emission-to-absorption efficiency for the gas opacity $\epsilon^{gas} = 1$. Dust will be hotter than gas in optically thin irradiated regions. In Section 8 we will briefly address again the effect of the gas opacity on the rim structure and the disk evolution.

### 3. THE STRUCTURE OF THE RIM

In this section we present the typical structure of the rim. We introduce model S100, which is characterized by a constant surface density of $\Sigma = 100$ g cm$^{-2}$. For this model we neglect the accretion stress and heating. The stellar parameters are the ones of a typical Herbig Ae class star (van den Ancker et al. 1998), with $T_{*} = 10,000$ K, $R_{*} = 2.5$ R$_{☉}$, and $M_{*} = 2.5$ M$_{☉}$. The resulting luminosity is $L_{*} = 56$ L$_{☉}$. We use a logarithmically increasing radial grid size, ranging from 0.2 to 4 AU with a total of 1280 grid cells. The vertical extent is 0.36 radians (20.6°) and is composed of 128 grid cells. The parameters of model S100 are summarized in Table 1.

The top panel of Figure 1 presents the convergence of model S100. The plot shows radial optical depth for the irradiation $\tau_{*} = 1$ at the midplane versus iteration number. The positions of three midplane temperatures, 1200, 800, and 400 K, are overplotted. Initially, there is no dust present and the disk is totally optically thin. In the first five iterations we increase

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top panel: convergence of the S100 radiation hydrostatic model. The vertical axis is the iteration number. At each iteration, the radius of unit starlight optical depth in the midplane is shown by a black vertical bar. Over the first five iterations, the amount of dust is increased logarithmically to its final value. Crosses, circles, and triangles mark where the midplane temperature is 1200, 800, and 400 K, respectively. Green vertical lines divide the inner disk into the dust-free zone (A), dust halo (B), inner rim (C), and shadowed region (D). Middle panel: final midplane radial temperature profile (black solid line) in model S100. Red curves show the temperatures of optically thin gas (dotted) and dust (dashed). The blue dotted line shows the dust sublimation temperature. Bottom panel: final midplane radial profile of the dust-to-gas mass ratio.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** 2D profiles of temperature (top) and dust density (bottom) in the R-Z/R plane, for the converged radiation hydrostatic model S100. The black lines indicate the optical depth unity for the irradiation (dashed line) and for the thermal emission (solid line). The yellow lines show the 1200 K (dashed) and 800 K (solid line) temperature contours. The green bars mark the positions of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$, which are defined by $\tau_{*} = 1$ at the midplane, and the position in which the incidence angle becomes zero.</figcaption></figure>

logarithmically the value of the dust-to-gas mass ratio to its final value. The results show that the radiation hydrostatic model quickly converges to a stable rim and temperature structure. We note that the optical depth from the inner gas disk in front of our computational domain (between 0.2 AU and three stellar radii) is $\tau_0 = 0.14$ at the midplane.

The final converged profiles of the temperature and dust-to-gas mass ratio are presented in Figures 1 and 2. The radial midplane temperature is plotted in the middle panel of Figure 1. Here we overplot the dust sublimation temperature from Equation (9) and the optically thin temperature of the dust and gas. The optically thin temperature for a given value of absorption-to-emission efficiency $\epsilon$ can be calculated with

$$T_{\text{thin}} = \left(\frac{1}{\epsilon}\right)^{0.25} \left(\frac{R_{*}}{2r}\right)^{0.5} T_{*}. \qquad (15)$$

The bottom panel of Figure 1 shows the midplane radial profile of the dust-to-gas mass ratio. The 2D profiles of the temperature and the dust density in the R-Z/R<fnref n="5" /> plane are plotted in Figure 2.

In the converged solution we define four distinct regions, marked with green vertical lines in Figure 1. Region A is the dust-free disk inward of 0.3 AU. In this zone the temperature follows the optically thin temperature of the gas ($\epsilon^{\text{gas}} = 1$). In zone B, the dust starts slowly to condense at the rate at which the temperature drops. This region can also be seen as an optically thin dust halo, a region with a small dust amount.

The inner dust halo is clearly visible in the bottom panel of Figure 2, as well as the rapid increase of the dust density at around 0.5 AU. We emphasize that this optically thin dust halo appears as a natural outcome of the radiation hydrostatic models. In this dust halo, the temperature is higher than the optically thin gas temperature due to the lower $\epsilon$ value of the dust; see Equation (15).

The border between regions B and C marks the region in which most of the dust suddenly condenses, which is the actual beginning of the rim. This is the position where the irradiation optical depth passes $\tau_{*} = 1$; see the top panel of Figure 1. The final radial position of $\tau_{*} = 1$, between zones B and C, depends on the global shape of the rim and hence the strength of backwarming, as well as the column of dust and gas in front of the rim. We call this point $R_{\text{rim}}^{\text{in}}$ from now on.

The rim occupies region C. Here the optical thickness quickly increases, while the temperature quickly drops. To understand the temperature profile across region C, we must consider the global shape of the rim, shown in Figure 2. The temperature is roughly constant between 0.5 and 0.8 AU. This is due to the high incidence angle at which the starlight strikes the rim surface. We define the rim's outer edge as the innermost point where starlight travels parallel to the surface of unit starlight optical depth. This point is marked in Figure 1 by the green vertical line dividing zones C and D, and in Figure 2 by a green bar. It corresponds to a local maximum in the aspect ratio Z/R of the infrared photosphere $\tau_{Z}^{\text{NIR}} = 1$. We call this point $R_{\text{rim}}^{\text{out}}$ from now on. We note that for model S100, the 1200 and 800 K contour lines are close to $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$. The locations of these midplane temperatures compared to the positions of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$ change on the order of 10% for the other models presented in this work.

<footnote n="5">⁵ The R-Z/R plane has the advantage that it shows the geometry of the rim along stellar rays (horizontal lines correspond to the path of individual photons from the star). The true geometry of the rim in the R-Z plane is illustrated in Figure 5 for the more realistic, uniform M models.</footnote>
<table><caption>**Table 2**
Model Input Parameter and Results</caption><thead><tr><th>Model</th><th>κ<sup>*</sup><sub>P</sub></th><th>R<sub>*</sub><br>R<sub>☉</sub></th><th>M<sub>*</sub><br>M<sub>☉</sub></th><th>T<sub>*</sub></th><th>M&#775;</th><th>L<sub>*</sub><br>L<sub>☉</sub></th><th>R<sup>in</sup><sub>rim</sub></th><th>R<sup>out</sup><sub>rim</sub></th><th>τ<sup>NIR</sup><sub>Z</sub> = 1/τ<sup>*</sup><sub>r</sub> = 1</th><th>h/r</th><th>P<sub>max</sub></th></tr></thead><tbody><tr><td>S100</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>-</td><td>56</td><td>0.46</td><td>0.92</td><td>0.10/0.14</td><td>0.04</td><td>...</td></tr><tr><td>MDe-8</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.42</td><td>0.89</td><td>0.12/0.15</td><td>0.04</td><td>0.85</td></tr><tr><td>MDe-7</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-7</sup></td><td>56</td><td>0.34</td><td>0.86</td><td>0.14/0.17</td><td>0.04</td><td>0.85</td></tr><tr><td>MDe-9</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-9</sup></td><td>56</td><td>0.57</td><td>0.93</td><td>0.10/0.14</td><td>0.04</td><td>0.87</td></tr><tr><td>LS21</td><td>1917</td><td>2.12</td><td>2.0</td><td>8500</td><td>10<sup>-8</sup></td><td>21</td><td>0.13</td><td>0.48</td><td>0.10/0.13</td><td>0.03</td><td>0.46</td></tr><tr><td>LS11</td><td>1784</td><td>2.0</td><td>1.8</td><td>7500</td><td>10<sup>-8</sup></td><td>11.3</td><td>0.09</td><td>0.31</td><td>0.09/0.11</td><td>0.02</td><td>0.28</td></tr><tr><td>α<sub>out</sub> = α<sub>in</sub></td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.35</td><td>0.77</td><td>0.09/0.13</td><td>0.04</td><td>...</td></tr><tr><td>α<sub>out</sub> = 10<sup>-4</sup></td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.43</td><td>0.93</td><td>0.14/0.17</td><td>0.04</td><td>0.93</td></tr><tr><td>T<sub>MRI</sub> = 800 K</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.41</td><td>0.97</td><td>0.12/0.15</td><td>0.04</td><td>0.97</td></tr><tr><td>T<sub>MRI</sub> = 1200 K</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.37</td><td>0.76</td><td>0.12/0.15</td><td>0.04</td><td>0.43</td></tr><tr><td>RHD_MDe-8</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-8</sup></td><td>56</td><td>0.43</td><td>0.91</td><td>0.12/0.15</td><td>0.04</td><td>0.90</td></tr><tr><td>RHD_MDe-9</td><td>2100</td><td>2.5</td><td>2.5</td><td>10000</td><td>10<sup>-9</sup></td><td>56</td><td>0.60</td><td>0.95</td><td>0.10/0.14</td><td>0.04</td><td>0.90</td></tr></tbody></table>

**Note.** Model name, irradiation opacity in cm$^2$g$^{-1}$, stellar radius, stellar mass, stellar surface temperature in Kelvin, mass accretion rate in units of solar mass per year, stellar luminosity, starting position of the rim position in units of AU, outer rim position in units of AU, maximum height of the NIR/optical surface at $R_{\text{rim}}^{\text{out}}$ in units of R, scale height $h/r = c_s/v_\phi$ at $R_{\text{rim}}^{\text{out}}$, and position of pressure maximum in units of AU.

Zone D is the region shadowed by the inner rim and starts near 1 AU in this model. Here the temperature drops below the optically thin gas temperature. We note again that the shadowing (zero incidence angle) can be seen in Figure 2, following the straight irradiation $\tau_* = 1$ line between 0.9 and 2 AU (green vertical bar). Radially outward of 2 AU, the disk starts to flare again (nonzero incidence angle).

In summary, the rim consists of three zones: (1) a hot, optically thin dust halo; (2) the starlit rim with its triangular cross section, whose radial extent is several times the density scale height; and (3) a shadowed, cool zone beyond the rim. All the models in this work display qualitatively similar structures despite widely differing parameters.

## 4. MODEL WITH A CONSTANT $\dot{M}$

We now move to the steady-inflow model with the surface density profile obtained as described in Section 2.3. In addition to the model described in the previous section, the gas surface density at each iteration is calculated according to Equation (12). We note that in the steady-inflow models, the accretion heating associated with the finite viscosity is neglected.

We first consider model MDe-8 with a typical mass accretion rate of $\dot{M} = 10^{-8} M_{\odot}$ per year. Previous works have shown that low or even vanishing accretion rates are able to fit the median SED of Herbig-type systems (Mulders & Dominik 2012). The results are presented in Table 2 and Figure 3. Both gas and dust surface densities show roughly power-law profiles with a jump at $T_{\text{MRI}}$. There, the surface densities increase by roughly one order of magnitude at the location where $\alpha$ drops by roughly one order of magnitude. In the annulus with increasing surface density, the gas is rotating at super-Keplerian velocities. The pressure maximum is located close to $R_{\text{rim}}^{\text{out}}$ (Figure 3, cross symbol). This is the location where we expect solid material to accumulate. The dust surface density strongly increases at the position of $R_{\text{rim}}^{\text{in}}$, while otherwise the dust surface density scales with the gas surface density. Overall, the rim structure looks very similar to that presented in Section 3. In model MDe-8, the rim's maximum angular extent as seen from the star is set by the annulus at 0.9 AU, where the starlight-absorbing surface rises to $Z/R = 0.15$ and the NIR-emitting surface to $Z/R = 0.12$.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Gas (solid line) and dust (dotted line) surface density radial profile for model MDe-8. The green vertical thick bars mark the positions of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$ to show the rim's extent. The cross marks the location of the pressure maximum.</figcaption></figure>

These values are similar to radiation hydrostatic models of the inner rim by Mulders & Dominik (2012) ($Z/R = 0.14$) or by Vinković (2014) ($Z/R = 0.11$). Table 3 presents the values of all the constants used in this paper.

## 5. SENSITIVITY TO THE FREE PARAMETERS

There are a number of free parameters that characterize the physical model. The most important ones are the value of the mass accretion rate and the stellar luminosity, which we will investigate in Sections 5.1 and 5.2, respectively. The effects of $\alpha_{\text{out}}$, the angular momentum transport rate in the dead zone, and the threshold temperature, $T_{\text{MRI}}$, are investigated in Sections 5.3 and 5.4, respectively.

### 5.1. The Influence of $\dot{M}$

To study the influence of the mass accretion rates, we add two models MDe-9 and MDe-7 with uniform accretion rates of $\dot{M} = 10^{-9}$ and $\dot{M} = 10^{-7} M_{\odot}$ per year. The results are summarized in Table 2 and Figure 4. Overall, the disk structure
<table><caption>**Table 3**
Physical Constants</caption><thead><tr><th>Constants Value</th><th>Name</th></tr></thead><tbody><tr><td>μ<sub>g</sub> = 2.353</td><td>Mean molecular weight</td></tr><tr><td>k<sub>B</sub> = 1.3806 × 10<sup>-16</sup> erg K<sup>-1</sup></td><td>Boltzmann constant</td></tr><tr><td>u = 1.6605 × 10<sup>-24</sup> g</td><td>Atomic mass unit</td></tr><tr><td>G = 6.6726 × 10<sup>-8</sup> cm<sup>3</sup> g<sup>-1</sup> s<sup>-2</sup></td><td>Gravitational constant</td></tr><tr><td>a<sub>R</sub> = 7.5657 × 10<sup>-15</sup> erg cm<sup>-3</sup> K<sup>-4</sup></td><td>Radiation constant</td></tr><tr><td>σ<sub>b</sub> = 5.6704 × 10<sup>-5</sup> erg cm<sup>-2</sup> s<sup>-1</sup> K<sup>-4</sup></td><td>Stefan-Boltzmann constant</td></tr><tr><td>c = 2.99792 × 10<sup>8</sup> cm s<sup>-1</sup></td><td>Speed of light</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Gas (solid line) and dust (dotted line) surface density radial profile for three values of the mass accretion rate $\dot{M}$. In each case, the dust surface density has been multiplied by a factor of 100. The green vertical thick bars mark the positions of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$ to show the rim's extent. The crosses mark the location of the pressure maximum. The green dotted line shows the critical surface density for which the disk becomes gravitationally unstable.</figcaption></figure>

at the rim looks similar for the three different mass accretion rates. Both gas and dust surface densities scale with $\dot{M}$. With increasing mass accretion rate, the position of $R_{\text{rim}}^{\text{in}}$ moves radially inward, mainly due to the higher optical depth $\tau_0$ and the higher vapor partial pressures in Equation (9). Special care should be taken when interpreting the results for the hydrostatic model $\text{MD}e-7$. Here we had to reduce the value of $\tau_0$ by a factor of 10 to prevent the inner disk ($r < 0.3$ AU) from becoming optically thick, which would move the rim toward small radii and so out of the computational domain. In addition, for this model the accretion heating is important (see Section 6).

Figure 4 includes the upper limit profile of the surface density at which the gravitational instability is triggered. We define the unstable surface density by

$$ \Sigma_Q = 1 = \frac{c_s^{\text{thin}} \Omega}{\pi G}, \qquad (16) $$

with the Toomre parameter $Q$ and the sound speed at optical thin dust temperature $c_s^{\text{thin}}$. As an example, a mass accretion rate of $\dot{M} = 10^{-7} M_{\odot} \text{yr}^{-1}$ and an accretion stress such that $\alpha_{\text{out}} = 10^{-4}$ would render the disk gravitationally unstable at about 5 AU (see Figure 4).

## 5.2. The Influence of the Stellar Luminosity

It is well known that the stellar luminosity is the most important parameter that determines the sublimation front location (Dullemond & Monnier 2010). We examine its effect on our results by calculating the structure of disks surrounding stars less luminous than considered above. More specifically, for the case $\dot{M} = 10^{-8} M_{\odot}$ per year, we add the following two models. The first uses $L_* = 21 L_{\odot}$ (for which we have $R_* = 2.12 R_{\odot}$, $M_* = 2.0 M_{\odot}$, and $T_* = 8500$ K), and the second, $L_* = 11.3 L_{\odot}$ (in which case $R_* = 2.0 R_{\odot}$, $M_* = 1.8 M_{\odot}$, and $T_* = 7500$ K). We named these models LS21 and LS11, respectively. Due to the different stellar surface temperature, we have to modify the dust opacities to $\kappa_p(T_*) = 1917 \text{ cm}^2 \text{ g}^{-1}$ and $1784 \text{ cm}^2 \text{ g}^{-1}$ and the computational domain radial extent to 0.1–2 AU and 0.07–1.5 AU, respectively. All other parameter values are kept fixed. The results are summarized in Table 2.

As expected, the rim disk structure clearly moves radially inward when the star luminosity decreases. For example, the location of $R_{\text{rim}}^{\text{in}}$ moves from 0.42 AU for a luminosity of $56 L_{\odot}$ to 0.09 AU for a luminosity of $11.3 L_{\odot}$ in model LS11. In addition, the decrease in stellar mass leads to an overall thinner disk at the position of $R_{\text{rim}}^{\text{out}}$. As this position, the scale height $h/r = c_s/v_{\phi}$ gradually decreases from 0.04 for the fiducial case to 0.033 and 0.026. Figure 5 summarizes the true geometry of the rim's surface in the $R-Z$ plane for the previous models.

## 5.3. The Influence of $\alpha_{\text{out}}$

In this section we investigate the effect of the value of $\alpha_{\text{out}}$ on the rim structure. To do so, we use the parameters of model $\text{MD}e-8$ and consider in addition the two cases given by $\alpha_{\text{out}} = \alpha_{\text{in}}$ and $\alpha_{\text{out}} = 10^{-4}$ (we remind the reader that $\alpha_{\text{out}} = 0.001$ in model $\text{MD}e-8$). The first case corresponds to a fully turbulent disk (i.e., without a dead zone), and the second describes a nearly laminar dead zone. The results are summarized in Table 2 and Figure 6. The model $\alpha_{\text{out}} = 10^{-4}$ yields results that are similar to model $\text{MD}e-8$. The rim radius $R_{\text{rim}}^{\text{in}}$ moves only slightly outward (from 0.42 to 0.43 AU). Changes remain modest but are more important for the case $\alpha_{\text{out}} = \alpha_{\text{in}}$. Here the most important difference is the absence of a pressure maximum. $R_{\text{rim}}^{\text{in}}$ decreases from 0.42 to 0.35 AU because of the reduced backwarming of the dust wall.

Overall, we conclude that the effect of varying $\alpha_{\text{out}}$ remains modest. This is because the rim is mainly located where $T > 1000$ K, whereas the surface density changes due to varying $\alpha_{\text{out}}$ occur at locations where $T < 1000$ K.

## 5.4. The Influence of $T_{\text{MRI}}$

For the last parameter comparison we investigate the effect of the value of $T_{\text{MRI}}$ on the rim structure. To do so, we use the parameters of model $\text{MD}e-8$ and consider in addition the two cases given by $T_{\text{MRI}} = 800$ K and $T_{\text{MRI}} = 1200$ K. The results are summarized in Table 2. The case $T_{\text{MRI}} = 800$ K is similar to model $\text{MD}e-8$. The outer rim radius $R_{\text{rim}}^{\text{out}}$ moves slightly outward (from 0.89 to 0.97 AU) due to the surface density increase at larger radii. The same happens with the position of the pressure maximum, which moves outward from $P_{\text{max}} = 0.85$ to 0.97 AU.

The case $T_{\text{MRI}} = 1200$ K has a stronger effect on the rim profile. In this case, the higher surface density at higher
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Shape of the irradiation optical depth unity surface for different mass accretion rates and stellar luminosities in the R–Z plane. We conserve the correct aspect ratio to show the true geometry. The vertical thick bars mark the position of $R_{\text{rim}}^{\text{out}}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Radial surface density profile of the gas (solid line) and dust multiplied by a factor of 100 (dotted line), for the models with modified $\alpha_{\text{out}}$. The vertical thick bars mark the positions of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$. The crosses mark the location of the pressure maximum.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

temperatures shifts the rim from $R_{\text{rim}}^{\text{in}} = 0.42$ to 0.37 AU and from $R_{\text{rim}}^{\text{out}} = 0.89$ to 0.76 AU. The pressure maximum shifts even more, from 0.85 AU in model MDe-8 to 0.43 AU. For model T<sub>MRI</sub> = 1200 K, we set the computational domain radial extent to 0.2–3 AU. In Section 7.2, we compare the previous models with the observational contraints.

## 6. RADIATION HYDRODYNAMICAL SIMULATIONS

Finally, we investigate the stability of the rim, using long-term 2D radiation hydrodynamical (RHD) simulations. The hydrostatic models are used as initial conditions. Our goals are twofold: in addition to verifying that the disk remains in steady state for thousands of dynamical times, we also want to investigate the effect of accretion heating on the disk structure. The numerical setup is almost identical to that used by Flock et al. (2013) and is briefly recalled in Appendix E, along with a more detailed description of the modifications required for the simulations to be completed.

We chose models MDe-8 and MDe-9 to match the mass accretion rate of typical Herbig star models (Mulders & Dominik 2012). The (RHD) simulations RHD_MD1e-9 and RHD_MD1e-8 were integrated for a runtime of 10,000 inner

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Top: time evolution of the radial surface density profile for the radiation hydrodynamical model RHD_MD1e-9. Bottom: radial surface density profiles of the radiation hydrostatic models MD1e-8 and MD1e-9 (solid lines) compared to the final snapshots of the radiation hydrodynamical models RHD_MD1e-8 and RHD_MD1e-9 (dashed lines). The vertical thick bars mark the position of $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$.</figcaption></figure>

orbits. In both simulations, the final disk structure after 10,000 inner orbits is almost identical to the initial state. In both models, the radial velocity fluctuations remain small, of the order of $10^{-3}c_s$ in the midplane. For illustrative purposes, we
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Radial midplane temperature profile of model RHD\_MD1e-8 (black solid line) and model RHD\_MD1e-8 (black dashed line). The red lines correspond to the optical thin temperature of the gas (red dotted line) and the dust (red dashed line). The blue line shows the sublimation temperature of the dust.</figcaption></figure>

present in Figure 7 (top panel) a spacetime diagram showing the evolution of the surface density radial profile for model RHD\_MD1e-9. The final surface density profiles of both models are compared in Figure 7 (bottom panel) to the initial static disk solution. The results show that the final profiles are very close to a passive disk solution (i.e., without accretion heating), especially for model RHD\_MD1e-9. For model RHD\_MD1e-8, the surface density adapts to a slightly different equilibrium as the effect of accretion heating becomes visible. The radial midplane temperature of model RHD\_MD1e-8 at this final state is plotted in Figure 8, overplotting also the initial temperature profile. The effect of the accretion heating becomes visible especially in the shadowed region. However, we note that overall the density and temperature structure of model RHD\_MD1e-8 remains similar to that of the radiation hydrostatic model MD1e-8.

## 7. COMPARISON WITH OBSERVATIONAL CONSTRAINTS

In this section we compare our results with observational constraints. We start by presenting the SEDs of our models in Section 7.1. We investigate the effect of the model parameters on the SED in Section 7.2. Finally, we compare our results for the rim radius with observational constraints in Section 7.3 and construct synthetic images in Section 7.4.

### 7.1. Spectral Energy Distribution

We use the Monte Carlo radiative transfer code RADMC3D (Dullemond 2012) to calculate the SED associated with the disk structures computed above. We consider seven individual wavelengths between 1 and 7 µm to cover the regime corresponding to our domain size and temperature range. For the dust opacity we use the wavelength-dependent table as shown in Figure 14 in Appendix A. For the calculation of the SED, we assume a disk inclination of 45° and scale the results to a distance of 122 pc (which corresponds to the distance of the Herbig star HD 163296). We compare our results with the Herbig star sample described by Mulders & Dominik (2012), for which model LS21 matches the Herbig A6 median stars' luminosity, surface temperature, and stellar mass (see

<figure><img src="image_12.png" /><figcaption>**Figure 9.** SED of the sample of Herbig stars (median spectral type A6) by Mulders & Dominik (2012) (black lines) with the upper and lower quartile (gray shade). The SED of model LS21 is shown with the pure dust component (solid blue line) and including an additional gas component (dashed line). The blue crosses show the seven individual wavelengths calculated for our models.</figcaption></figure>

Appendix B therein). We calculate the SED using two strategies. In the first, we consider the contribution of the dust component only. In the second, we add the gas component, assuming a gray gas opacity. The results of both models are plotted in Figure 9 and compared with the sample from observations. The NIR emission of model LS21 is below the observational median, especially at 2 µm, where the emission is lower by a factor of two. The gas only contributes a few percent of additional flux in this regime.

### 7.2. Effect of Accretion Heating, $\alpha_{out}$, and $T_{MRI}$ on the SED

In this section we investigate the effects of the accretion heating, $\alpha_{out}$, and $T_{MRI}$ on the SED. Overall the effects remain small. We compare the SED of the hydrostatic model MD1e-8 with the profile of the radiation hydrodynamical model RHD\_MD1e-8, using the simulation output we obtained at the last time step. The two resulting SEDs are very close, with differences less than 5%.

Similarly small changes are seen for model $\alpha_{out} = 10^{-4}$, which has a higher dust density at the outer rim position. For the case $\alpha_{out} = 10^{-4}$ the emission at 4 µm is greater by 8% than the model MD1e-8 with $\alpha_{out} = 0.001$. For model $\alpha_{out} = \alpha_{in}$, there is no surface density increase at the outer rim and the emission at 4 µm decreases by 17% compared to the model with $\alpha_{out} = 0.001$. Comparing with the observational constraints, the models would favor a strong drop of $\alpha_{out}$, but the changes remain small and we cannot rule out the possibility of a constant $\alpha_{out}$ from the SED alone.

Finally, model $T_{MRI} = 800$ K shows a 6% decrease of emission at 4 µm as the surface density jump moves radially outward, affecting the emission at longer wavelengths. Model $T_{MRI} = 1200$ K shows a small increase of 5% at 2 µm but a larger decrease of 10% at 5 µm compared to model MD1e-8 with $T_{MRI} = 1000$ K.

### 7.3. Rim Radius

As we have seen, the sublimation front takes a pointed shape spread over several pressure scale heights in radius. It is thus
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Inner dust ring radius over luminosity relation. The squares, circles, and rectangles represent T Tauri, Herbig Ae, and Herbig Be stars, respectively, adapted from Dullemond & Monnier (2010) and Millan-Gabet et al. (2007, p. 539). We overplot our results by colored vertical lines, showing the radial extent $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$ in Table 2. The different models are LS11 (red), LS21 (green), MDe-9 (yellow), MDe-8 (blue), and MDe-7 (magenta). We note that for the MDe-X models, the value of $R_{\text{rim}}^{\text{out}}$ is similar.</figcaption></figure>

not straightforward to define the exact radial position of the rim. In this section, we compare the observational determined radii presented by Millan-Gabet et al. (2007, p. 539) and Dullemond & Monnier (2010) with the rim radial extent $R_{\text{rim}}^{\text{in}}$ and $R_{\text{rim}}^{\text{out}}$ from our models (see Table 2). As shown in Figure 10, the theoretical and observational rim radii compare very well for the range of stellar luminosity we explored. This comparison highlights once more that the actual rim becomes more radially extended with higher mass accretion rate, as shown by the vertical length of the bars in Figure 10 for the MDe-X models. By contrast, the actual ring position is determined by the central star luminosity and agrees with the systematic trend seen in the observations. We note that we did not include the effect of accretion luminosity, as this remains small for the considered systems (Mulders & Dominik 2012).

Overall, the SED and the position of the rim from our models are in broad agreement with previously published models by Dullemond & Monnier (2010) and Mulders & Dominik (2012), which do not include the gas between the rim and the star. As previous models have shown, there is still a lack of emission at 2 µm wavelength by roughly a factor of two. Neither the presence of a small dust halo, nor the contribution of the gas, nor accretion heating, nor changing the surface density or the position of the surface density increase is able to solve that problem. As recently shown by Turner et al. (2014), the effect of magnetic pressure, namely, a thickening of the disk at the rim location, remains a possible solution.

## 7.4. Synthetic Images

To finally provide a realistic view of the rim, we constructed synthetic images of our models. Figure 11 shows the central region of model RHD_MD1e-8 in steady state for an inclination of 45°. The three panels from top to bottom are synthetic images at wavelengths 1.25, 2.2, and 4.8 µm (J, K, and M bands). The field of view is 2 AU wide and includes the hot optical thin gas and the actual rim (compare zones B and C in Section 3). The J and K bands both show a bright inner ring at

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Synthetic images of the final output from the radiation hydrodynamical model RHD_MD1e-8, viewed 45° from face-on. The intensity maps correspond to 1.25 (top), 2.2 (middle), and 4.8 µm (bottom).</figcaption></figure>

the $\tau_{*} = 1$ midplane region, while the extended triangular cross section of the rim is cooler and emits more in M band. The small dust halo in front of the rim emits slightly in the J band; see Figure 11, top panel. Figure 12 shows the three bands combined into a color image, with J, K, and M mapped to the blue, green, and red channels, respectively. The plots shows that overall the K and M emission is strongest (red-yellow colors), while the contribution from J bands remains small.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Synthetic color image of the final output from the radiation hydrodynamical model RHD_MD1e-8, viewed 45° from face-on. The blue, green, and red channels come from the 1.25, 2.2, and 4.8 µm bands, respectively. The three channels share a common, linear intensity scale determined by the extrema of the 4.8 µm image.</figcaption></figure>

## 8. DISCUSSION

This work represents a first step in constructing self-consistent models of the inner regions of protoplanetary disks that account for both dynamical and thermodynamical constraints while remaining numerically tractable. There are still some important limitations to this work, which we will overview in the following.

We considered a uniform dust-to-gas mass ratio of 0.01. However, turbulence, dust settling, and radial drift all change the local dust density and consequently the opacity. In addition, small grains could be quickly depleted due to the fast growth and settling (Brauer et al. 2008; Birnstiel et al. 2010; Zsom et al. 2011; Okuzumi et al. 2012). Our models show a pressure maximum appearing at the location of the ionization temperature $T_{\text{MRI}}$ at which the surface density increases due to the drop of accretion stress. At that location, larger dust particles could be concentrated and increase the dust-to-gas mass ratio. Collisions between the larger particles could provide the small grains, which could be mixed in the upper layers. It is possible that such an increase in dust density could help increase the height of the rim and create a larger shadowed region. However, due to the fact that this location is at lower temperatures ($T_{\text{MRI}} < T_{\text{ev}}$), it is located further outward in the disk and is unlikely to increase the flux at 2 µm. In addition, the rim can also be affected by the photoelectric heating (Thi et al. 2011) and the radiation pressure on the dust (Vinković 2014). A more sophisticated treatment of the dust is needed in the future.

In this work, we have fixed the gas opacity to a small value of $\kappa_{\text{gas}} = 10^{-4} \text{ cm}^2 \text{ g}^{-1}$. For this value, the inner gas disk remains optically thin for a given range of mass accretion rates from $\dot{M} = 10^{-9}$ to $\dot{M} = 10^{-8} M_{\odot}$ per year. One way to change the gas opacity without affecting the rim position would be to change the value of the accretion stress in the ionized region.

Assuming an accretion stress of $\alpha_{\text{in}} = 0.1$ for temperatures above 1000 K would reduce the gas surface density by one order of magnitude and thus allow a higher gas opacity of $\kappa_{\text{gas}} = 10^{-3}$ without affecting the optical depth of the inner gas disk. Finally, we note again that a detailed implementation of gas line radiation transfer and the frequency-dependent gas opacity would go far beyond the scope of this work.

## 9. CONCLUSIONS

We have developed the first radiation hydrodynamical models of the silicate sublimation front in protoplanetary disks around Herbig Ae stars. The models are axisymmetric and include stellar irradiation, dust and gas opacity, dust sublimation, and condensation. The effects of turbulence (angular momentum transport and dissipative heating) are modeled by means of a temperature-dependent kinematic viscosity. This dependence is chosen to capture the onset of magnetorotational turbulence due to dust grains' thermionic emission at temperatures above about 1000 K (Desch & Turner 2015). The models are inflow-equilibrium solutions with radially constant mass accretion rates. We compute cases with mass flow rates from $\dot{M} = 10^{-9}$ to $10^{-7} M_{\odot}$ per year and stellar luminosities from $L_{*} = 11 L_{\odot}$ to $L_{*} = 56 L_{\odot}$.

Using numerical simulations in which we solve the time-dependent radiative-hydrodynamics equations, including viscosity, we have shown the stability of the rim to axisymmetric modes. In addition, the models remain in steady state for thousands of dynamical timescales. For $\dot{M} \le 10^{-8} M_{\odot}$ per year, accretion heating has only a modest effect on the flow. The height of the NIR optical depth unity $\tau_{Z}^{\text{NIR}} = 1$ is at around $Z/R = 0.1$ to $Z/R = 0.14$, depending on the model parameters and the mass accretion rate. The SEDs of our models, calculated by Monte Carlo radiative transfer tools, show good agreement with previous radiation hydrostatic models (Mulders & Dominik 2012; Vinković 2014).

The inner disk can be divided into three structures, summarized in Figure 13. The first is an optically thin halo of hot dust between the sublimation front and the star. Such hot dust in front of the rim has been proposed in previous models (Kama et al. 2009), though its extent and optical depth are lower than predicted by Vinković et al. (2006). The second structure is the actual rim front, where the dust condenses over the geometrically thin layer, absorbing most of the starlight. This front has a triangular cross section in the poloidal plane, with the point lying in the midplane aimed at the star. The triangular cross section is quite radially extended, stretching several times the gas density scale height. The third structure is the material lying outside the sublimation front, which is shadowed from the starlight by the optically thick dust. The radial position of the rim in our models matches well with observational constraints from NIR interferometry. These results depend on the parameters as follows:

1. The radial extent of the inner rim triangular cross-section shape depends mainly on the mass accretion rate. This extent could be several tenths of an AU. Higher mass accretion rates increase this radial extent of the rim. Such an extent is consistent with NIR interferometry observations, which have clearly shown a smooth, radial extended structure, derived from visibility curves (Tan-nirkulam et al. 2008; Benisty et al. 2010). Finally, even though the mass accretion rate has an effect on the radial
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Schematic structure of the inner rim of the disk around a young intermediate-mass star, summarizing the results of our radiation hydrodynamical models. This meridional cross section has the star and the system's rotation axis at left. Vertical dotted lines divide the disk into, from left to right, (1) optically thin gas with silicate vapor; (2) the silicate sublimation front, forming the inner rim of the optically thick, dusty material—the starlight absorbed at the front heats the rim, making it vertically extended; (3) the shadow cast by the rim; and (4) the outer disk, whose slightly upward-curving surface lets it too see the star. The local pressure maximum just beyond the rim is a location where planet-forming solids can accumulate at temperatures near 1000 K.</figcaption></figure>

extent, it has only a small effect on the SED and the actual radial position of the rim.

2. We find that the accretion heating only affects the temperature profile in the shadowed region behind the inner rim for the given parameter space we investigated. We observe no substantial increase of the actual inner rim height with increasing accretion rate. For our models, accretion heating becomes important for a mass accretion rate above $\dot{M} \ge 10^{-8} M_{\oplus}$ per year. The relative weak dependence on accretion heating can also be explained by the fact that the high accretion stress is expected for temperatures above the ionization temperature of 1000 K and such regions are optically thinner than the regions below 1000 K.

3. For the given ionization temperature we found that the location of the pressure maxima is at a region with temperatures around 1000 K. Solids are likely to concentrate near this radius under gas drag forces, since the pressure gradients make the gas rotate faster than Keplerian inside the maximum and slower than Keplerian outside (Haghighipour & Boss 2003). Concentration of pebbles or boulders followed by collisional fragmentation could lead to more abundant submicron dust, which would increase the height where the optical depth falls to unity. At this region, the temperatures are sufficiently high for annealing to form crystalline silicates.

In summary, we have presented the first models that enable us to study the inner dust rim region with radiation hydrodynamical simulations. The work should be seen as a bridge, connecting previous highly sophisticated dust Monte Carlo radiation transfer hydrostatic models of the inner rim and global hydrodynamical simulations of stratified disks. With a simplified radiation transfer, we capture the relevant physics,

and the presented models compare very well with the observational constraints of ring radius and SED. The NIR emission at 2 µm is about half as bright relative to the host star as observed, similar to existing models (Mulders et al. 2010). Whether the missing NIR flux can come from a magnetically supported atmosphere remains an open question (Turner et al. 2014). Better characterizing the environments near young stars of different masses is essential if we are to understand the origins of the population of close planets.

For the first time we have performed axisymmetric radiation hydrodynamical simulations of this region, and we confirm the flow's stability over thousands of orbits. Furthermore, robust dust abundance and temperature distributions can be obtained on grids coarse enough that 3D calculations are now feasible.

The authors thank Antonella Natta, Wlad Lyra, Rafael Millan-Gabet, Gijs Mulders, and Satoshi Okuzumi for useful comments on the manuscript. We thank Andrea Mignone for supporting and advising us with the newest PLUTO code. Parallel computations have been performed on the Genci supercomputer "curie" at the calculation center of CEA TGCC and on the zodiac supercomputer at JPL. For this work, S.F. and M.F. received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement nr. 258729. This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration and with the support of the NASA Exoplanet Research program via grant 14XRP14_20153. Copyright 2016 California Institute of Technology. Government sponsorship acknowledged.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Wavelength dependence of the dust opacity. The dust opacity by Draine & Lee (1984) (blue) and Preibisch et al. (1993) (red) is plotted for comparison. The two dust opacities used for model MDe-8 are overplotted.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Converged radial temperature profiles for different sublimation functions. The red curves correspond to the optical thin temperature of the gas (red dotted line) and the dust (red dashed line). The blue dotted line shows the sublimation temperature of the dust.</figcaption></figure>

# APPENDIX A
## DUST OPACITY

The dust opacities are calculated by Mie theory using the method by Wolf & Voshchinnikov (2004). To calculate the opacity, we assume a range of silicate and carbon particles between $a_{min} = 5$ µm and $a_{max} = 100$ µm with a size distribution profile $\sim a^{-3.5}$ and a silicate abundance of 62.5%. The profile of the dust opacity is plotted in Figure 14. For comparison, we plot also opacity values by Preibisch et al. (1993) and by Draine & Lee (1984), which assume slightly smaller dust sizes. To reduce the complexity of the problem, we use gray opacities for the simulations. For example, we define the Planck mean opacity as

$$ \kappa_P(T) = \frac{\int \kappa_{\nu} B(\nu, T) d\nu}{\int B(\nu, T) d\nu}, \qquad (17) $$

with the Planck function $B(\nu, T)$. We note that in this work we assume the same evaporation temperature for silicate and carbon grains. Especially refractory carbon grains could survive to higher temperatures. We will address this in a future work.

# APPENDIX B
## TEST OF THE DUST SUBLIMATION FUNCTION
### FOR $T > T_{ev}$

We perform additional models in which we replace the function for $T > T_{ev}$ in Equation (10) with an exponential function as by Kama et al. (2009),

$$ f^{\text{k09}} = \exp\left\{-a_{\text{steep}}\left(\frac{T_{\text{ev}} - T}{T}\right)^2\right\}, \qquad (18) $$

with the steepness factor $a_{steep}$. In the following, we show the influence of the functions in Equation (10) on the converged dust temperature profile. The first function for $T > T_{ev}$ causes a fast sublimation of the dust for temperatures higher than the sublimation temperature. This function has only a minor influence on the final temperature profiles, compared to different functions with different steepness; see Figure 15.

Independent of our choice of evaporation function, the models converge always to a solution with an optical thin dust halo in front of the rim (compare the region with higher temperature in front of the rim). The function for $T < T_{ev}$ in Equation (10) smooths out the high opacity jump due to the dust condensation and so allows us to resolve the irradiation absorption at the rim front independent of the grid resolution and optical thickness. An example of a calculation without using this function is presented in a later section. We also note that the value of $f_{\Delta\tau}$ effects only the steepness of the temperature drop between the $\tau_{*} = 1$ and $\tau_{*} = 100$ positions. In Figure 15 we perform additional radiation hydrostatic versions of model S100, following Equation (18) and using two different steepness factors of $a_{steep} = 400$ and $a_{steep} = 1600$. The plot shows that the functions have only a minor effect on the location of the rim. In addition, independent of the steepness factor, we naturally obtain an optically thin dust halo in front of the rim, which can be seen at the region with high temperature $T > T_{thin}^{\epsilon=1}$.

# APPENDIX C
## TEST OF THE DUST SUBLIMATION FUNCTION
### FOR $T < T_{ev}$

In this test, we want to show the importance of the dust sublimation function Equation (10) for $T < T_{ev}$. We performed the model RHD_S100, which uses as initial conditions the results of model S100. We also performed a model without the second function in Equation (10), here called model RHD_S100<sup>MOD</sup>. Figure 16 compares the two snapshots of the 2D temperature profile of model RHD_S100 and model RHD_S100<sup>MOD</sup> after 200 inner orbits. In model RHD_S100<sup>MOD</sup>, the irradiation heating is absorbed in one cell, which causes jumps in the temperature profile. We note that this effect becomes even stronger for simulations with higher surface density and so higher optical depth. With this test we want to underline the importance of resolving the irradiation absorption and therefore the use of Equation (10).
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Comparison of 2D temperature profile in the R-Z/R plane, after six inner orbits for model RHD_S100 (left) and model RHD_S100<sup>MOD</sup> (right). The green lines indicate the optical depth unity for the irradiation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Temperature profiles of model MD1e-8 using full frequency-dependent irradiation and temperature-dependent dust opacity (dashed line) compared to the model using a constant two-opacity model (solid line) over radius (left) and height (right).</figcaption></figure>

## APPENDIX D
### FREQUENCY-DEPENDENT IRRADIATION AND TEMPERATURE-DEPENDENT DUST OPACITY

In this appendix we compare our simplified constant two-
opacity model, described in Section 2 with a more complex
model, using frequency-dependent irradiation and Planck mean
dust opacity, which depends on temperature $\kappa_P(T)$. For further
details, e.g., on the calculation of the frequency-dependent
irradiation, we refer to our previous work (Flock et al. 2013).
For this additional model we use 60 frequency bins to sample
the opacity, presented in Figure 14. In addition, the Planck and
Rosseland opacities are calculated in each cell for the given
dust amount and dust temperature. The resulting midplane and
vertical temperature profiles for the two models are presented in
Figure 17. There is no significant difference visible between
those two models, in terms of both temperature and structure of
the rim. The fixed opacity model reproduces very well the
structure and temperature profile of the more complex model
using frequency-dependent irradiation. We note that in the
fixed opacity model, we overestimate the local Planck opacity
in a large area in the disk by a factor between 1 and 2.

## APPENDIX E
### MODIFICATIONS OF THE DYNAMICAL RHD SETUP

The 2D radiation hydrodynamical simulations we present in
Section 6 were performed using a second-order in space and

time numerical configuration of the PLUTO code (Mignone
et al. 2012). We solve the same RHD equations as in Flock
et al. (2013), with two modifications: the magnetic field is set to
zero, and we include in the momentum and total energy
equations the effect of a finite kinematic viscosity. Its
amplitude is calculated as described in Section 2.3. We used
the Harten–Lax–Van Leer Riemann solver with a Courant
number of 0.3 to increase the numerical stability.

The domain size is set to $R_{\text{in}}-R_{\text{out}} = 0.3-3.0$ au and $\Delta\theta = 0.36$ radians, and we used $1024 \times 128$ grid cells, respectively, in the radial and meridional directions. The initial conditions are provided by the radiation hydrostatic models MDe-9 and MDe-8 and are interpolated on the computational grid using the built-in interpolation routine in PLUTO, which allows us to change the domain size and grid resolution.

Three modifications were necessary to increase the numer-
ical stability of the simulations:

1. A buffer zone extending from 0.3 to 0.32 AU and from 2.5 to 3.0 AU in radius, where we reestablish the surface density on a timescale shorter than the dynamical timescale. This prevents the loss of material due to the pure outflow radial boundary condition.

2. A time adaptation routine for the dust sublimation, which smooths out strong fluctuations in the local dust amount for each time integration step. As a result, changes of the dust density are limited to 10% per time step. This
<figure><img src="image_21.png" /><figcaption>Figure 18. Radial midplane temperature (top left) and radial dust-to-gas mass ratio (top right) for different resolutions. Bottom: 2D temperature contour plot for the highest-resolution model (10240 × 1024 cells).</figcaption></figure>

prevents sudden changes in the opacity that could lead to problems for the radiative solver to converge.

3. A modified gravitational potential in the radial inner and vertical upper layers. This was done by modifying the gravitational potential according to

$$ \Phi = \frac{GM_*}{r - 0.175(0.4 \text{ au} - r)(\theta - \pi/2)^2} \quad \text{for } r < 0.4 \text{ AU} \qquad (19) $$

$$ \Phi = \frac{GM_*}{r} \quad \text{elsewhere.} \qquad (20) $$

Its effect is to reduce the density contrast from 13 to 9 orders of magnitude. We found that this small change in the uppermost layers substantially increases the stability of the dynamical calculations.

## APPENDIX F
RESOLUTION TEST

To check whether the rim structure is robust against changes in the spatial resolution, we perform a convergence study on model S100. The three new hydrostatic solutions have resolutions 2, 4, and 8 times finer than the model described in Section 3. Figure 18 shows that the higher the resolution, the thinner the layer near unit starlight optical depth, where the dust abundance is limited by Equation (10), and the steeper the

temperature gradient. However, over this eightfold range in spatial resolution, the temperatures away from the front vary by at most a few percent. We conclude that despite the smoothing applied to the dust abundance near the rim, our calculations recover the correct overall structure.

## REFERENCES

* Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1

* Barge, P., & Sommeria, J. 1995, A&A, 295, L1

* Benisty, M., Natta, A., Isella, A., et al. 2010, A&A, 511, A74

* Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tuscon, AZ: Univ. of Arizona Press), 691

* Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79

* Bitsch, B., Morbidelli, A., Lega, E., Kretke, K., & Crida, A. 2014, A&A, 570, A75

* Boley, A. C., Morris, M. A., & Ford, E. B. 2014, ApJL, 792, L27

* Bowler, B. P., Johnson, J. A., Marcy, G. W., et al. 2010, ApJ, 709, 396

* Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

* Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53

* Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., et al. 2001, ApJ, 547, 1077

* Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52

* Decampli, W. M., Cameron, A. G. W., Bodenheimer, P., & Black, D. C. 1978, ApJ, 223, 854

* Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156

* Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

* Dullemond, C. P. 2012, RADMC-3D: A Multi-Purpose Radiative Transfer Tool, Astrophysics Source Code Library, ascl:1202.015

* Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

* Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205
Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70

Faure, J., Fromang, S., Latter, H., & Meheut, H. 2014, arXiv:1411.3236

Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43

Flock, M., Henning, T., & Klahr, H. 2012, ApJ, 761, 95

Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343

Haghighipour, N., & Boss, A. P. 2003, ApJ, 598, 1301

Helling, C., Winters, J. M., & Sedlmayr, E. 2000, A&A, 358, 651

Hillenbrand, L. A., Strom, S. E., Vrba, F. J., & Keene, J. 1992, ApJ, 397, 613

Hu, X., Zhu, Z., Tan, J. C., & Chatterjee, S. 2015, arXiv:1508.02791

Ilee, J. D., Fairlamb, J., Oudmaijer, R. D., et al. 2014, MNRAS, 445, 3723

Isella, A., & Natta, A. 2005, A&A, 438, 899

Isella, A., Testi, L., & Natta, A. 2006, A&A, 451, 951

Johnson, J. A., Morton, T. D., & Wright, J. T. 2013, ApJ, 763, 53

Kama, M., Min, M., & Dominik, C. 2009, A&A, 506, 1199

Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21

Kraus, S. 2015, Ap&SS, 357, 97

Kretke, K. A., & Lin, D. N. C. 2012, ApJ, 755, 74

Kretke, K. A., Lin, D. N. C., Garaud, P., & Turner, N. J. 2009, ApJ, 690, 407

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

Lissauer, J. J., Dawson, R. I., & Tremaine, S. 2014, Natur, 513, 336

Lyra, W. 2014, ApJ, 789, 77

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 479, 883

Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009, A&A, 497, 869

Lyra, W., & Mac Low, M.-M. 2012, arXiv:1204.5711

Lyra, W., Paardekooper, S.-J., & Mac Low, M.-M. 2010, ApJL, 715, L68

Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91

Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478

Matsumura, S., Pudritz, R. E., & Thommes, E. W. 2009, ApJ, 691, 1764

McClure, M. K., D'Alessio, P., Calvet, N., et al. 2013, ApJ, 775, 114

Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476

Mendigutía, I., de Wit, W. J., Oudmaijer, R. D., et al. 2015, MNRAS, 453, 2126

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Millan-Gabet, R., Malbet, F., Akeson, R., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tuscon, AZ: Univ. of Arizona Press), 539

Millan-Gabet, R., Schloerb, F. P., & Traub, W. A. 2001, ApJ, 546, 358

Mulders, G. D., & Dominik, C. 2012, A&A, 539, A9

Mulders, G. D., Dominik, C., & Min, M. 2010, A&A, 512, A11

Muzerolle, J., D'Alessio, P., Calvet, N., & Hartmann, L. 2004, ApJ, 617, 406

Natta, A., Prusti, T., Neri, R., et al. 2001, A&A, 371, 186

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106

Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615

Preibisch, T., Ossenkopf, V., Yorke, H. W., & Henning, T. 1993, A&A, 279, 577

Reffert, S., Bergmann, C., Quirrenbach, A., Trifonov, T., & Künstler, A. 2015, A&A, 574, A116

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94

Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117

Tannirkulam, A., Monnier, J. D., Harries, T. J., et al. 2008, ApJ, 689, 513

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tuscon, AZ: Univ. of Arizona Press), 339

Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711

Turner, N. J., Benisty, M., Dullemond, C. P., & Hirose, S. 2014, ApJ, 780, 42

Umebayashi, T., & Nakano, T. 1988, PThPS, 96, 151

van den Ancker, M. E., de Winter, D., & Tjin A Djie, H. R. E 1998, A&A, 330, 145

Varnière, P., & Tagger, M. 2006, A&A, 446, L13

Vinković, D. 2012, MNRAS, 420, 1541

Vinković, D. 2014, A&A, 566, A117

Vinković, D., Ivezić, Ž, Jurkić, T., & Elitzur, M. 2006, ApJ, 636, 348

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Wolf, S., & Voshchinnikov, N. V. 2004, CoPhC, 162, 113

Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73