Planet formation and migration near the silicate sublimation front in protoplanetary disks

The increasing number of newly detected exoplanets at short orbital periods raises questions about their formation and migration histories. A particular puzzle that requires explanation arises from one of the key results of the Kepler mission, namely the increase in the planetary occurrence rate with orbital period up to 10 days for F, G, K and M stars. We investigate the conditions for planet formation and migration near the dust sublimation front in protostellar disks around young Sun-like stars. For this analysis we use iterative 2D radiation hydrostatic disk models which include irradiation by the star, and dust sublimation and deposition depending on the local temperature and vapor pressure. We perform a parameter study by varying the magnetized turbulence onset temperature, the accretion stress, the dust mass fraction, and the mass accretion rate. Our models feature a gas-only inner disk, a silicate sublimation front and dust rim starting at around 0.08 au, an ionization transition zone with a corresponding density jump, and a pressure maximum which acts as a pebble trap at around 0.12 au. Migration torque maps show Earth- and super-Earth-mass planets halt in our model disks at orbital periods ranging from 10 to 22 days. Such periods are in good agreement with both the inferred location of the innermost planets in multiplanetary systems, and the break in planet occurrence rates from the Kepler sample at 10 days. In particular, models with small grains depleted produce a trap located at a 10-day orbital period, while models with a higher abundance of small grains present a trap at around a 17-day orbital period. The snow line lies at 1.6 au, near where the occurrence rate of the giant planets peaks. We conclude that the dust sublimation zone is crucial for forming close-in planets, especially when considering tightly packed super-Earth systems.


Introduction
Our knowledge of the origins of planetary systems relies on our understanding of the disks of gas and dust orbiting young stars. With the growing number of detected exoplanetary systems, we want to learn about their formation and evolution, starting from their building blocks in young protoplanetary disks. One key result of the Kepler mission was the discovery of many planetary systems containing Earth-sized and super-Earth planets orbiting with periods between one day and one year (Lissauer et al. 2011;Batalha et al. 2013;Fabrycky et al. 2014). The Kepler mission focused on searching for planetary systems around F-, G-, and K-type stars. Further analysis of the results showed a clear rise in the occurrence rate of super-Earths (rocky planets of several Earth masses) with orbital periods longer than 10 days. The occurrence is uniform at longer periods (Youdin 2011;Mulders et al. 2015;Petigura et al. 2018). Over recent years, possible mechanisms for this rise were investigated, including truncation of the disk by the stellar magnetosphere (Lee & Chiang 2017), trapping of planets at the inner edge of the disk (Brasser et al. 2018;Carrera et al. 2019), and the removal of close-in planets by tidal interactions inside 10-day orbital periods (Rice et al. 2012). Even more recently, Mulders et al. (2018) have found, when including observational constraints on multiplanetary systems, that on average the innermost planet orbits with a 12-day period.
Broadly speaking, models of the formation of short-period planets are based either on the drift of pebbles to locations near the inner edge of the disk where they concentrate and grow into planets (Boley et al. 2014;Chatterjee & Tan 2014) or on planets that form farther out in the disk, perhaps by pebble accretion (Lambrechts et al. 2019;Izidoro et al. 2019), and then migrate inward before reaching locations near the inner edge of the disk where migration halts because the disk-planet torques go to zero (Terquem & Papaloizou 2007;Ida & Lin 2010;Faure & Nelson 2016;Izidoro et al. 2017). To compare these planet formation models with observations, the main disk property used is the abundance of solids (Dawson et al. 2015(Dawson et al. , 2016. Formation models are typically based on simple assumptions about the inner disk structure. More detailed self-consistent disk profiles of density and temperature including irradiation and dust sublimation have not yet been computed.
In this work our aim is to develop irradiated hydrostatic disk models around young Sun-like stars. With respect to the Kepler mission results, we want to address the influence of the disk structure on both pebble drift and planet migration, and in particular show how the temperature and density structure leads to locations where drifting pebbles and migrating planets are trapped.
We follow the numerical methods that we used previously to model the inner disk around Herbig Ae/Be-type stars (Flock et al. 2016). In this work, we have adapted all the parameters to match a protoplanetary disk around a young Sun-like star, which is less massive and less luminous than Herbig Ae/Be stars, and for which the disks sustain a higher mass accretion rate.
We first investigate the density and thermal structure of the disk, and its impact on pebble drift and planet migration. We consider steady-state models with a radially uniform mass accretion rate. To determine the surface density we use the temperature dependent α value, which we directly calculate from our previous 3D radiation magneto-hydrodynamic (MHD) simulations (Flock et al. 2017). In our previous work we found that the effect of viscous heating is small (Flock et al. 2016), and therefore we focus on passively heated disks in this work.
For temperatures below 850-1000 K (which is dependent on the density), the ionization level drops (Desch & Turner 2015) and turbulent activity due to the magneto-rotational instability (MRI) quickly decreases. Among the main processes governing the temperatures in the region close to the star are the irradiation heating and sublimation of silicate grains, which change the opacity by orders of magnitude (Pollack et al. 1994), and which we include in these models.
The structure of the paper is as follows. In Sect. 2 we briefly review the main parts of the numerical method. In Sect. 3 we present the 2D radiation hydrostatic solutions. In Sect. 4 we calculate migration maps. Sections 5 and 6 follow with the discussion and conclusions.

Method
We followed the numerical method presented by Flock et al. (2016). We first calculated a radiation hydrostatic equilibrium, solving iteratively the radiation transfer including irradiation. In the first step the gas surface density profile was determined using at the cylindrical radius R, assuming a constant radial mass accretion rateṀ. The viscosity is given by where c s is the sound speed, Ω = GM * /R 3 is the disk rotation frequency, and the stress-to-pressure ratio is given by with α MRI being the stress-to-pressure ratio in the active zone (T > T MRI ), α DZ being the stress-to-pressure ratio in the dead zone, and T MRI being the ionization transition temperature for the MRI to operate. The initial temperature field T (r, θ) was calculated using the optically thin solution, which depends on the stellar luminosity and the opacities. We then calculated the density ρ(r, θ) and azimuthal velocity v φ (r, θ) by solving the governing equations in spherical coordinates (r, θ, φ) with axisymmetry to obtain the solution for hydrostatic equilibrium, where Φ is the gravitational potential Φ = GM * /r (with G being the gravitational constant and M * the stellar mass) and P is the thermal pressure that relates to the temperature through the ideal-gas equation of state: with the mean molecular weight µ g , the Boltzmann constant k B , and the atomic mass unit u. For a given density field, the internal energy and radiation field in radiative equilibrium were obtained as the steady-state solution to the following coupled pair of equations describing heating, cooling, and the flux-limited diffusion of the disk's thermal radiation: Here we used the adiabatic index γ, the radiation energy E R , the irradiation flux F * , the flux limiter λ (Levermore & Pomraning 1981), the Rosseland and Planck mean opacity κ R and κ P , the radiation constant a R = 4σ b /c (with the Stefan-Boltzmann constant σ b and the speed of light c). We computed F * by solving the transfer equation along radial rays from the central star, treated as a point source. Following our previous work, we simplified the problem by assuming κ R = κ P and considering three opacities: the gas opacity κ gas , the Planck dust opacity at the dust sublimation temperature κ P (T s ), and the Planck dust opacity at the stellar temperature κ P (T * ).
We took the gas opacities from Malygin et al. (2014) for starlight and for the disk's thermal emission. The two dust opacities were frequency-averaged, one over the blackbody spectrum at the sublimation temperature and the other over the irradiating starlight spectrum. The combined opacity was determined using the dust-to-gas ratio. For details on the gas and dust opacities we refer to the appendix.
We ignored the heat released by accretion within the disk. Viscosity was only used to set the constant mass accretion and to obtain a surface density profile consistent with the uniform mass accretion rate. In our previous work we found the effect of viscous heating remains small (Flock et al. 2016), while recent disk wind models predict an even smaller viscous heating effect (Béthune et al. 2017;Mori et al. 2019) in the dead-zone region. We address this again in the discussion section.
Convergence was reached by iteratively solving for the radiation equilibrium and hydrostatic equilibrium, including irradiation from the star. We iterated until we reached convergence. We iterated 30 times (compared to 20 times in our previous method) as these models span a larger domain.

Stellar parameters
T * = 4300 K, R * = 2.6 R M * = 1.0 M Opacity κ P (T * ) = 1300 cm 2 g −1 κ P (T s ) = 700 cm 2 g −1 κ gas ∼ 10 −5 cm 2 g −1 Grid size N r × N θ 2048 × 128 Cell aspect ratio ∆R/R∆θ ∼ 1.1 Radial domain R in −R out 0.05-10 AU Vertical domain θ π/2 ± 0.15 Ionization transition T MRI = 900 K α viscosity α MRI = 10 −2 for T > T MRI α DZ = 10 −3 for T ≤ T MRI Dust-to-gas mass ratio f 0 = 10 −3 of small grains (≤10 µm) Mass accretion rateṀ = 3.6 × 10 −9 M yr −1 As in our previous work, we followed Pollack et al. (1994) and used the fitting model of Isella & Natta (2005) that applies to situations for which the most refractory grains are silicates when determining the sublimation temperature: The value of T s was then used to calculate the dust-to-gas ratio f D2G , which is the ratio between the dust density ρ d and the gas density. We slightly modified our previous formula (Eq. (11) in Flock et al. 2016) to with the additional constraint that f D2G ≤ f 0 , with f 0 being the maximum dust-to-gas mass ratio and f ∆τ = 0.2/(ρ dust ∆rκ P (T * ) being the dust amount to account for optical depth of ∆τ = 0.2 for a given grid cell with radial width of ∆r. This value ensures that we resolve the absorption of the incoming radiation at the inner rim. We slightly steepened the transition compared to our previous work where we used ∆τ = 0.3. We also centered the absorption at τ = 2/3 where most of the absorption takes place.
For T < T s and τ > 3.0 most of the irradiation is absorbed and the available silicates occur in solid form. For more information on the method, the flux limiter, and the treatment of the optical depth at the rim wall we refer to Flock et al. (2016).

Model parameters
The model input parameters are summarized in Tables 1 and 2.
The stellar parameters were taken from stellar evolutionary tracks by Siess et al. (2002) andD'Antona &Mazzitelli (1994), and represent a young solar-mass star with an age of approximately one million years. Here we recalculated the stellar spectrum-weighted irradiation opacity. The redder spectrum of Table 2. Model name, ionization transition in Kelvin, MRI α value, α in the dead zone, mass accretion rate in (M yr −1 ), dust-to-gas mass ratio, gas opacity in (cm 2 g −1 ), and resulting location of the pebble trap in au.
The dust-to-gas mass ratio ranges from 10 −4 to 10 −3 representing grains with sizes up to 10 µm, which dominate the opacity at the relevant temperatures and wavelengths. Modeling indicates that most of the solid material is locked in larger grains (Birnstiel et al. 2012), consistent with a protostellar disk millimeter flux-radius correlation (Rosotti et al. 2019) and the weakness of silicate features in the infrared spectra of the disks around young stars with masses near solar (Furlan et al. 2006(Furlan et al. , 2009(Furlan et al. , 2011. We investigated ionization transition temperatures of 900 and 1000 K. For our reference model we assume 900 K, motivated by the ionization thresholds, which are at cooler temperatures for higher densities (Desch & Turner 2015). We studied accretion rates between 10 −9 and 10 −8 M yr −1 corresponding to those of disks in the Chamaeleon I star-forming region (Manara et al. 2017). For our reference model, MREF, we used a value of 3.6 × 10 −9 M yr −1 . We note that as we did not include the outer edge of the disk in our models, the total disk mass remains a free parameter.
For the stress-to-pressure ratio α MRI in the MRI active region, we took values between 0.01 and 0.1. The larger value was determined from our previous 3D radiation non-ideal MHD simulations. In Flock et al. (2017) we found a stress-to-pressure ratio of up to 0.1 for the net vertical magnetic flux case. For temperatures below the ionization threshold we chose stress-topressure ratios α DZ between 5 × 10 −4 and 10 −3 , to mimic the accretion activity either by hydrodynamical instabilities (Lyra & Umurhan 2019) or a magnetically driven wind (Béthune et al. 2017). We note again that the increase in the stress-to-pressure ratio at the ionization transition has a direct impact on the surface density profile. Compared to our previous models (Flock et al. 2016) we also investigated larger jumps in α and the surface density.
The grid resolution is 2048 cells in radius logarithmically spaced and 128 in the meridional direction for a domain extending from 0.05 to 10 AU in radius and to elevation angles 0.15 radians on either side of the equatorial plane in the meridional direction. The inner radial boundary was chosen to be close to but still outside the radius where the stellar magnetic field is expected to truncate the disk. A147, page 3 of 13 A&A 630, A147 (2019) Fig. 1. Temperature in the meridional plane for our reference model MREF, a young solar-type star with a radially-uniform mass accretion rate of 3.6 × 10 −9 M yr −1 , a dust-to-gas mass ratio for small grains of 0.001 and stress-to-pressure ratios of α MRI = 0.01 and α DZ = 0.001. The green solid line shows optical depth τ = 2/3 for the starlight irradiation. The yellow dashed line shows the silicate sublimation front. The white solid line shows the water sublimation front. The white circle shows the location of the pebble trap, where large grains would concentrate. Fig. 2. Density in the meridional plane of model MREF, a young solar-type star with a radially uniform mass accretion rate of 3.6 × 10 −9 M /yr, a dust-to-gas mass ratio for small grains of 0.001 and stress-to-pressure ratios of α MRI = 0.1 and α DZ = 0.001. Contours of the gas and dust density are shown for 10 −9 and 10 −12 g/cm 3 (yellow solid line). The green solid line shows the optical depth line τ = 2/3 for the starlight irradiation. The white solid line shows the water snow line in the disk. The white circle shows the location of the pebble trap.

Results
The results derived from the radiation hydrostatic solution for our reference model are shown in Figs. 1 and 2. Figure 1 presents the 2D temperature profile of the disk, showing a similar structure to that found for a Herbig-type star (Flock et al. 2016), but where the key features occur on smaller radial scales closer to the star: a hot dust halo in front of the rim inside 0.08 au, a curved dust rim between 0.08 and 0.15 au, a small shadowed region between 0.2 and 0.3 au, and a flared disk beyond 0.3 au. The starlight optical depth τ = 2/3 line starts at the midplane at 0.08 au, while reaching an elevation angle of θ = 0.08 rad at 0.15 au. At 8 au the τ = 2/3 line reaches θ = 0.15, which corresponds to roughly 3 scale heights. We find a pressure maximum and pebble trap located at 0.12 au, indicated by a white circle in Fig. 1  drifting pebbles are trapped, and could possibly trigger planetesimal formation. Figure 2 presents the corresponding 2D gas density profile. The next feature with increasing distance from the star is the rise in gas density at 0.13 au due to the lower stress-to-pressure ratio at temperatures below the ionization threshold. Midplane densities of 10 −8 gcm −3 are reached in the gas. The density jump is clear in the gas surface density profile in the top panel of Fig. 3. For this model we obtain a gas surface density of 680 g cm 2 at 1 au. Figures 3  and 4 show the radial profile of the disk scale-height-to-radius ratio (H/R) and the midplane temperature. Near the silicate sublimation front the H/R is around 0.02, generally increasing outward to reach 0.028 at 1 au. The disk radial temperature shows the typical profile, with a shallower decrease at the rim at around 1000 K, while at 1 au the temperature reaches roughly 200 K. For comparison, Figs. 3 and 4 also show the slopes of r −1 for density, r 1/4 for H/R, and r −0.5 for the temperature.
Figures 1 and 2 both include the snow line. We locate the snow line using the water vapor mole fraction In Fig. 1 we determine the position where half of the water is in the form of ice, such that x ice = 0.5, using Using this formula we determine that the snow line crosses the midplane at 1.6 au, reaching 2 au at 2 scale heights above the midplane. Due to the lower pressure and higher temperature, the snow line moves radially outward with increasing vertical height until it reaches 7 au at 3 scale heights. The snow line is an important location for planet formation. In particular, due to the predicted agglomeration of larger grains this is a preferred location of giant planet formation (Oka et al. 2011;Ros & Johansen 2013;Drążkowska & Alibert 2017;Schoonenberg et al. 2018). We note that for our choice of stellar parameters and resulting stellar luminosity of ∼1.92 L the snow line is located farther out radially compared to models with a lower stellar luminosity. This value is similar to recent solar-mass and metallicity stellar models (Amard et al. 2019), which predict roughly 2 L for a star one million years in age.

Migration maps
In this section we calculate the migration torques on planets embedded in our reference model. We are especially interested in determining where orbital migration halts, as this may set up the architecture of the new planetary system.
We use the torque formula from Paardekooper et al. (2011) as adapted by Bitsch & Kley (2011). The total torque Γ tot is the sum of the Lindblad torque Γ L and the corotation torque Γ c . The Lindblad torque arises because of the spiral waves launched at the planet's Lindblad resonances. We follow Paardekooper & Papaloizou (2008) and write with α s the negative slope of the surface density profile Σ ∝ r −α s and β s the negative slope of the temperature profile T ∝ r −β s . The normalization torque Γ 0 is defined as with q the planet-to-star mass ratio, h the aspect ratio h = H/R, Σ P the surface density at the planet's location r P , and Ω P the orbital frequency of the planet. Calculating the corotation torque Γ c is more complicated because this can saturate at a level depending on planet mass, disk viscosity, and local thermal diffusion rate. We compute Γ c as detailed in Sects. 5.6 and 5.7 of Paardekooper et al. (2011). The above calculation of Γ tot is valid until the gap-opening or type II migration regime is reached. To determine the gapopening mass we follow the equation by Crida et al. (2006): with the Hill radius and the Reynolds number Planets exceeding this mass would substantially change the surface density structure, and hence the migration torque.
In Fig. 5 we present a map of the total torque in the planet orbital radius-mass plane. For this calculation, the values of the density, temperature, stress-to-pressure ratio, surface density, and opacity are all taken from the reference radiation hydrostatic solution which is in inflow equilibrium. The figure shows that across most of the inner disk and for most planetary masses, migration is inward, as expected from type I migration models (Ward 1997;Fogg & Nelson 2007;Crida & Bitsch 2017). From left to right, the first region of outward migration is connected to the dust sublimation zone at around 0.08 au and the resulting drop in temperature. This region is relevant only for planets more massive than the Earth and it could be a halting point for planets all the way up to the gap-opening mass.
The next and most prominent region of outward migration lies in the steep rise of the surface density starting at 0.1 au. Here outward migration occurs for planets with masses from as low as 1% of an Earth mass and up to 100 Earth masses. The normalized migration torque is very strong, reaching Γ tot /Γ 0 ≥ 15. Such a planet trap at surface density transition was already proposed by Masset et al. (2006). The third outward migration region is due to the small shadowed zone where the temperature drops significantly starting at 0.13 au. This outward migration region affects planets with masses between 0.1 and 4 Earths. The radii where the torque vanishes (black solid contour lines) are especially interesting for planet formation and evolution models. At these locations, the Lindblad and corotation torques exactly cancel each other. A147, page 6 of 13 In Fig. 5 we overplot with a vertical green dashed line the location of the pressure maximum where inward-drifting pebbles accumulate. This location approximately coincides with the planet migration trap. If a steady radial drift of pebbles is present, we might expect a first planetary core to form at this location. Growth by accreting pebbles could continue until the planet reaches the pebble isolation mass. Lambrechts et al. (2014) found the pebble isolation mass to be Recently, Bitsch et al. (2018) found that the pebble isolation mass can, for certain parameters, be a factor of two higher than the fit by Lambrechts et al. (2014). In Fig. 5 we show both estimates. From this calculation we expect the planets to form at the pebble trap near 0.13 au and to have masses between 1.5 and 3 Earth masses.
Finally, gray shading indicates the region of type II migration where the planet's tides would be strong enough to open a gap. Such a drastic change in the disk structure means the migration torque calculation is invalid in this region. We also note that we do not see any change in the migration zones radially outward between 1 and 10 au. In Sect. 4.3 we compare these results with the exoplanet occurrence rate obtained from the Kepler sample.

Parameter study
In this section we vary the parameters (see Table 2). A summary of the resulting temperature and surface density profiles is in Fig. 6. Increasing the value of α MRI in the inner disk in model M0 compared to MREF leads to a steeper gradient in the surface density. This also affects the planet trap, extending the zero-torque region to larger planet masses (see Appendix A). The position of the pebble trap is shifted outward by only 0.01 au. Decreasing the abundance of small grains, and hence the opacity, has a larger effect. Decreasing the dust-to-gas ratio to 10 −4 in model M1 shifts the inner rim closer to the star, and results in a steeper rise in the temperature at the rim. The pebble trap moves radially inward to 0.093 au, corresponding to a 10-day orbital period. The other noticeably different model is M4 with the highest abundance of small grains, which shifts the pebble trap to 0.15 au, corresponding to a 22-day orbital period. Increasing the gas opacity by one order of magnitude in model M7 had only a small effect on the structure as the gas between the sublimation front and the star remains optically thin.

Comparison with exoplanet occurrence rates
In this section we compare our results with previous analyses of exoplanet occurrence rates. We note that for all models we use a fixed stellar type, while the Kepler sample is based on a mix of F-, G-, and K-type stars.
In the top panel of Fig. 7 we show the inferred real distribution of the orbital period of the innermost planet in Kepler multiplanetary systems reported by Mulders et al. (2018) 1 . Overplotted are the positions of the pebble traps we found in our models.
The pebble trap in our models reproduces quite well the position of the innermost planet in multiplanetary systems. The pebble trap occurs slightly farther from the star, with a mean orbital period of around 17 days. The dust-depleted model, M1, matches the 10-day peak in the measured planets' distribution. Mulders et al. (2018) show that the innermost planets had between one and four Earth radii, roughly corresponding to between one and ten Earth masses. Interestingly, this is similar to the mass range of the type I migration trap (Fig. 5). The scatter in the measured radii (and masses, assuming an Earth-like density) of the innermost planets is consistent with the range from the pebble isolation mass up to the gap-opening threshold.
The innermost planet trap in Fig. 5 is located near an 8-day orbital period and connected with the dust sublimation front. Inward-migrating super-Earths can also pile up at the planet trap at the inner edge of the gas disk where a cavity forms because of the stellar magnetosphere. That inner-edge trap sits closer to the star than the trap in our models, which is a result of the transition between the inner MRI-active region and the dead zone. We further note that multiplanetary systems in resonant chains can push the innermost planets farther in Carrera et al. (2018Carrera et al. ( , 2019 and Izidoro et al. (2019).
In the bottom panel of Fig. 7 we show the results of Fernandes et al. (2019) for the occurrence rate of giant planets from RV surveys and Kepler. They find a clear peak around 2 au, which is very close to the predicted position of the snow line in our models.

Dependence on luminosity
The position of the sublimation front is determined by the stellar luminosity (Dullemond & Monnier 2010). The position of the pebble trap has a similar dependence, as shown in Fig. 8 where we combine model MREF with the results for higher luminosity stars from Flock et al. (2016). The location of the pressure maximum varies with stellar luminosity to the power 0.6.
As with the silicate sublimation front, the snow line is determined by the stellar luminosity. Bitsch et al. (2015) show that the snow line can lie inside 1 au at accretion rates similar to those we consider. The snow line is significantly farther from the star in our models for two reasons. First, the stellar luminosity is about twice as great. Second, the opacity for starlight exceeds that for the disk's own thermal emission, where Bitsch et al. (2015) had the reverse. Absorbing efficiently and re-emitting inefficiently makes the dust hotter at a given starlight flux, pushing the snow line outward.

Discussion
In this section we discuss the effects of viscous heating, the magnetically-driven wind, and the stellar magnetic field on our model. We compare our results to previously found migration maps and other planet formation models. Finally, we report on the effect that eccentricity and resonant chains might have.

Viscously heated or passive disk
We have focused here on irradiated passive disks. Our similar radiation hydrodynamical models of disks around intermediatemass young stars show little impact from accretion heating (Flock et al. 2016). However for T Tauri disks accretion is more important since the ratio of accretion flux to irradiation flux at the sublimation front varies roughly inversely with stellar luminosity. To examine whether accretion heating is important, in Fig. 9 we compare the passive irradiated disk with one heated only by accretion. Following Hubeny (1990) the accretion flux radiated from the disk surface is using the corotation radius R 0 . With the accretion power deposited in the interior, the midplane temperature is for vertical Rosseland mean optical depth τ Z . The midplane accretion temperature for the reference model's surface density profile is the blue dashed line in Fig. 9. The accretion heating shifts the edge of the dead zone from 0.1 au to 0.26 au. However, in the passive model (solid line in Fig. 9) the annulus between these two radii is cool enough to remain MRI inactive. Furthermore, the turbulent surface layer or base of the magnetized wind, where the accretion power is deposited when the interior is inactive, has an optical depth comparable to or less than unity (Bai & Goodman 2009;Béthune et al. 2017;Mori et al. 2019). According to Eq. (20) the accretion temperature at the base of the active layer is well below the irradiation temperature, so accretion heating is negligible. Thus, the accretion flow in the passive irradiated model does not significantly alter its temperature. Furthermore, for the purpose of Fig. 9 we ignore A147, page 8 of 13 the dust sublimation; by speeding up cooling sublimation would act as a thermostat and limit the temperature to T s , again bringing the temperatures close to those in the passive model. We also note that the accretion-heated disk is much colder than the passive disk inside 0.08 au, due to the low gas opacity. If both irradiation and accretion were included, the temperature here would be similar to the irradiation-only case.
In summary, the passive irradiated models accurately represent the thermal structure. However, accretion heating is important between the inner rim and the edge of the dead-zone for higher mass accretion ratesṀ ≥ 10 −8 , pushing the edge of the dead zone and the planetary trap farther out. However, even in this case, we still expect that the first planetary core to form will be dry, since the dead-zone pebble trap always occurs at high temperatures of around 800-1000 K.

Torque map comparison
One noticeable difference from previous models is the contribution of the entropy-related corotation torque. In our models this effect is slight, giving rise only to a small outward migration region that sits outside the pebble trap in Fig. 5. In previous studies for T Tauri systems that include a higher viscous heating (Coleman & Nelson 2014;Bitsch et al. 2015), the outward migration zone extends to 5 au and beyond because of the steeper temperature gradient obtained in viscously heated disks. Following the previous discussion, we expect that the contribution from the entropy-related corotation torque in the dead zone should remain small in the absence of another source of heating.

Gap-opening mass
Equation (15) yields traditional gap-opening mass appropriate to viscous disks. However, if α DZ takes the form of an external torque applied by magneto-centrifugally launching a wind from the disk surface, its character may be quite unlike a viscosity, influencing the gap-opening mass. In this case, where the disk midplane is laminar, the gap-opening mass is much lower (Rafikov 2002;Dong et al. 2011;McNally et al. 2019), especially for the low values of H/R in our models.

Exact shape and position of the inner rim
When the sublimation front is highly unresolved, the starlight is absorbed in the near face of the first dust-containing cell it meets, yet the heat is deposited uniformly through that cell, including its unlit interior and back side. Therefore the cell's front side is cooler than it should be. Yet the front side is what we see in the infrared if the cell is optically thick. Thus, the disk's thermal radiation appears at too long a wavelength and the front's shape and position may be computed incorrectly. To resolve this issue, we developed a method that reduces the dust abundance smoothly across the sublimation front so that the transition from optically thin to thick is spatially resolved (Flock et al. 2016). This approach may become unnecessary in future when multiple grain compositions and sizes are included. A mix of species, each with its own sublimation temperature, naturally blurs the rim (Kama et al. 2009).
The rim's location suffers from uncertainty because anything near the star that is opaque enough allows dust to survive within its shadow. Our models, which assume a very low gas opacity, should therefore be seen as producing an outermost limit with respect to the radial position of the rim. We note that the exact location of the pebble and planetary traps are influenced by the rim location and the location of the ionization transition. Furthermore, we expect the small grains in the dead zone to deplete quickly. Ueda et al. (2019) show that small grains grow efficiently in the dead zone, due to its weak turbulence, forming pebbles that drift inward quickly.
One set of parameters we have not varied is the stellar temperature, radius, and mass. These parameters affect the position of the inner rim and other features in the structure of the disk, and thus should be investigated in the future (Mulders et al. 2015).

Influence of stellar magnetic field on the rim shape
The star's magnetic field could potentially disturb the disk near the sublimation front (Romanova et al. 2012). To determine whether our neglect of the stellar field is valid, we consider a dipole with strength of 1 kilogauss at the stellar surface (Johns-Krull 2007). The field strength in the midplane at 0.1 au is 1.8 gauss. For our reference model MREF, the thermal pressure at the midplane in the same location is P th = 240 dyn cm −2 . The plasma beta is thus 1900. If it penetrated the disk's plasma, a field of this strength would be important for driving MRI turbulence. Although we would not expect a very strong effect on the rim profile or height, which we investigated in our previous work (Flock et al. 2017), we note that the α MRI = 0.1 value used in model M2 might be more realistic under these conditions.

Comparison with other planet formation models
Because our disk models produce both a pebble trap and a planet trap, our results in principle support both an in situ planet formation model based on the drift and concentration of pebbles which then form planets (Boley et al. 2014;Chatterjee & Tan 2014) and a model based on the formation of planets farther from the star, which then migrate inward until they are halted at the planet trap (Izidoro et al. 2017. The migration model appears most able to explain the fact that some systems, such as Trappist I which contains seven terrestial planets around an M star (Gillon et al. 2016(Gillon et al. , 2017, are composed of planets in chains of mean motion resonances. In this model the planets migrate toward the inner edge of the protoplanetary disk until they are parked in resonant chains. Recent work suggests that many of these resonant A147, page 9 of 13 A&A 630, A147 (2019) chains can become unstable after the gas disk dissipates, explaining why most of the systems discovered by the Kepler mission are not in resonance (Izidoro et al. 2017Lambrechts et al. 2019).
Without an inner disk edge and planet trap, migration would drive all planets into the central star, and resonant chains would not exist because planets could not pile up. Ormel et al. (2017) proposed that planetary cores form at the snow line and then migrate inward while continuing to grow via pebble accretion. It is noteworthy that the innermost Trappist-1 planet is located at around 0.01 au, which is roughly consistent with the inner edge of the dust disk if the luminosity was at least one order of magnitude higher than it is currently, as is expected for young, very low-mass stars (Laughlin & Bodenheimer 1993).

Eccentricity effect on the torque
We assume an orbital eccentricity of zero when calculating the torque map. However, planets undergoing type I migration may have a small eccentricity value, for example due to being embedded in a turbulent region of the disk (Laughlin et al. 2004) or because multiple planets migrate together with a small separation. In such a case the outward migration zone would shrink since the corotation torque becomes weaker (Fendyke & Nelson 2014;Coleman & Nelson 2016). This could shift the planet traps even closer to the star, depending on the details of the density profile near the ionization transition zone. The shift is in the right direction to improve the fit to the occurrence rate of Kepler planets. In addition, if the turbulence is particularly vigorous, the outward migration zone located around 0.08 au may shrink, or disappear altogether if the eccentricity exceeds the disk aspect ratio of about 0.02. This will need to be examined in the future using radiation-MHD simulations of planets embedded in magnetized disks.

Effect of resonant chains
The innermost planet in a multiplanetary system could also be pushed nearer to the star by torques from the outer planets in the resonant chain (Carrera et al. 2018(Carrera et al. , 2019Izidoro et al. 2019). This might improve the agreement between the innermost planet position and the Kepler results of Mulders et al. (2018). Furthermore, multiplanetary systems tend to become unstable, leading to scattering and collisions among the planets once the disk has dispersed, which could explain the occurrence of super-Earths interior to any planet trap position (Izidoro et al. 2017.

Conclusion
We developed radiation hydrostatic models to describe the thermal and density structure of protoplanetary disks around T Tauri stars. The models are 2D and axisymmetric, and include stellar irradiation, dust and gas opacities, dust sublimation, and condensation. Magnetically driven accretion is modeled by means of a temperature-dependent kinematic viscosity. This dependence is chosen to capture the onset of magneto-rotational turbulence at temperatures around 900 K. The models are inflow-equilibrium solutions with radially constant mass accretion rates. All our models are for a typical young Sun-like star with parameters of T * = 4300 K, R * = 2.6 R , and M * = 1.0 M . We examine disk mass accretion rates fromṀ = 10 −9 to 10 −8 M yr −1 , and dust-to-gas mass ratios between 10 −4 and 10 −3 , considering grains smaller than 10 µm that are responsible for absorbing incoming stellar radiation and reemitting it in the thermal infrared.
The computational domain spans the dust sublimation front and the water snow line. The resulting disk structure is approximately a scaled-down version of the Herbig disk models obtained by Flock et al. (2016): -In our reference model the dust sublimation front occurs at around 0.08 au in the midplane, and curves out to around 0.15 au at an elevation angle of Z/R = 0.08. -The next feature outward after the sublimation front is a steep increase in density due to the drop in the ionization level and the corresponding drop in the accretion stress. Here at 0.13 au we find a robust local pressure maximum capable of trapping pebble-sized particles. -Orbital migration torques calculated for embedded super-Earths indicate that type I migration halts at around 0.13 au.
Outward torques on super-Earths nearer the star are strong enough to hold back a resonant chain with two or more outer planets of similar mass. The planet migration trap lies close to the peak (from the Kepler survey) in the orbital distribution of the innermost planets in multiplanetary systems. -The orbital period at the pebble trap in our reference model is 17 days. In a model with the dust abundance ten times lower, the pebble trap is shifted inward to an orbital period of 10 days. -For all of our models, the snow line is located close to 1.6 au at the midplane, reaching 2 au at two scale heights above the midplane. This closely matches the peak in the occurrence rate of giant planets from radial velocity surveys. We point out that our results can in principle support both of the leading planet formation scenarios. The disk structure we determined allows inward-drifting pebbles to accumulate and form planets at the pressure maximum a short distance outside the dust sublimation front. The model also allows planet formation to occur farther out in the disk, for example at the snow line, after which the planets migrate inward and become trapped beneath the inner rim formed by the dust sublimation front. In the future it will be important to test these findings using full 3D radiation-MHD simulations with embedded pebbles and planets, so the dynamics of these bodies can be computed within the context of the complex and turbulent flows present in the inner regions of protoplanetary disks.
We show torque maps for all remaining models in Fig. A.1. They are directly comparable with the map for the reference model in Fig. 5. In our previous work (Flock et al. 2016) we included a simple temperature-and pressure-independent gas opacity of κ gas = 10 −4 cm 2 g −1 . For this new work, we make use of the more detailed gas opacities derived by Malygin et al. (2014). Figure B.1 shows the Rosseland mean gas opacities.
Over the temperature and pressure regime of our model we derive a mean value of about κ gas = 10 −5 cm 2 g −1 . For simplicity we assume the same gas opacity also for the irradiation. Future calculations of the detailed gas opacity including photochemistry are necessary to improve this step. A higher gas opacity could affect the rim shape and the temperature profile of the inner disk.

B.2. Dust opacity
To compute the dust opacity we used the tool MieX (Wolf & Voshchinnikov 2004). We assumed silicate and carbon particles with radii between a min = 0.005 µm and a max = 10 µm and a grain size distribution having exponent −3.5. The particles are 62.5% astrophysical silicate and 37.5% graphite. The wavelength dependence of the dust opacity is plotted in Fig. B.1.
From the wavelength-dependent opacity we determine the Planck mean using with the Planck function B(ν, T ). The frequency integral corresponds to a wavelength between 0.05 µm and 2 mm. For the starlight color temperature T * = 4300 K, we determine two significant figures κ P (T * ) = 1300 cm 2 g −1 . For the sublimation front color temperature T s = 1500 K, we determine κ P (T s ) = 700 cm 2 g −1 .