Planet-disk interaction at the dead zone inner boundary

Inside-Out Planet Formation. III. Planet-disk interaction at the dead zone inner boundary


The Kepler mission has discovered more than 4000 exoplanet candidates. Many of them are in systems with tightly packed inner planets (STIPs). Inside-Out Planet Formation (IOPF) (Chatterjee & Tan, 2014) has been proposed as a scenario to explain these systems. It involves sequential in situ planet formation at the local pressure maximum of a retreating dead zone inner boundary (DZIB). Pebbles accumulate at this pressure trap, which builds up a pebble ring, and then a planet. The planet is expected to grow in mass until it opens a gap, which helps to both truncate pebble accretion and also induce DZIB retreat that sets the location of formation of the next planet. This simple scenario may be modified if the planet undergoes significant migration from its formation location. Thus planet-disk interactions play a crucial role in the IOPF scenario. Here we present numerical simulations that first assess the degree of migration for planets of various masses that are forming at the DZIB of an active accretion disk, where the effective viscosity is undergoing a rapid increase in the radially inward direction. We find that torques exerted on the planet by the disk tend to trap the planet at a location very close to the initial pressure maximum where it formed. We then study gap opening by these planets to assess at what mass a significant gap is created. Finally we present a simple model for DZIB retreat due to penetration of X-rays from the star to the disk midplane. Overall, these simulations help to quantify both the mass scale of first, “Vulcan,” planet formation and the orbital separation to the location of second planet formation.

Subject headings:
protoplanetary disks, planet-disk interactions, planets and satellites: formation


1. Introduction

Since launch in 2009, Kepler has revealed more than 4000 exoplanet candidates (e.g., Mullally et al., 2015). A large percentage () of these candidates are in systems with tightly-packed inner planets (STIPs). These are systems that usually have 3 or more detected planets of radii and with periods less than 100 days (Fang & Margot, 2012). There are two main scenarios that can produce such close-in planets: (1) formation in the outer disk followed by inward migration (e.g., Kley & Nelson, 2012; Cossou et al., 2013, 2014); (2) formation in situ (Hansen & Murray, 2012, 2013; Chiang & Laughlin, 2013; Chatterjee & Tan, 2014, hereafter CT14).

The inward migration scenario tends to produce planets that are trapped in orbits of low order mean motion resonances, which is not a particular feature of STIPs (Baruteau et al., 2014; Fabrycky et al., 2014). Thus it has been proposed that the lack of resonant pile-ups might be explained by a lower efficiency of resonance trapping or breaking of resonance by later dynamical processes for the typically low-mass Kepler-detected planets (Goldreich & Schlichting, 2014; Chatterjee & Ford, 2015).

The in situ formation scenario faces the challenge of concentrating enough solids in the inner region (Raymond & Cossou, 2014; Schlichting, 2014). For example, supply of pebbles by radial drift may be truncated if planet formation, perhaps initiated by the pebble streaming instability (Youdin & Goodman, 2005) , occurs in the outer disk (Lambrechts & Johansen, 2014; Levison et al., 2015; Bitsch et al., 2015). In situ formation models also face the challenge of reproducing the observed mass versus orbital radius distributions once effects of gas on protoplanet migration during the oligarchic growth phase are accounted for (Ogihara et al., 2015).

The Inside-Out Planet Formation (IOPF) scenario proposed by CT14 is a new type of in situ formation model. It starts with pebble delivery to the midplane transition region between the innermost MRI-active zone and a nonactive “dead zone,” where there is a local pressure maximum. The pebbles trapped by the pressure maximum will build up in a ring, which then forms a protoplanet, perhaps involving a variety of processes including streaming (Youdin & Goodman, 2005), gravitational (Toomre, 1964) and/or Rossby wave (Varnière & Tagger, 2006) instabilities. The protoplanet is expected to continue its growth, especially by pebble accretion, until it becomes massive enough to open a gap in the disk. This gap pushes the pressure maximum outwards by a few Hill radii thus creating a new pebble trap that is displaced from the planet’s orbit (Lambrechts et al., 2014), but may also allow MRI-activation in the region beyond the planet, which could induce further outward retreat of the DZIB.

As discussed by CT14, planetary migration, either before gap opening (Type I) or after gap opening (Type II), may in principle alter this scenario. However, the expectation is that the pressure maximum associated with this inner disk transition region will also act as a planet trap (Masset et al., 2006; Matsumura et al., 2009; Kretke & Lin, 2012; Bitsch et al., 2014), so the planet will stay at its initial formation location and keep accreting more materials, at least until the point of gap opening. However, this needs to be confirmed for the particular disk structure involved in IOPF and this is one of the goals of this paper.

The disk studied by Masset et al. (2006) was a purely passive disk, i.e, the temperature was set only by irradiation from the central star. Zhang et al. (2014) studied planetary migration in a disk with an inner region heated by viscous dissipation and an outer region heated by stellar luminosity. This structure also contained a local pressure maximum that acted as a trap for Type I planetary migration. In the model for IOPF, the first planet is formed at the outer edge of the MRI active zone around the protostar (where  K, leading to thermal ionization of alkali metals) and migration needs to be studied in the context of an active disk, i.e., heated by viscous accretion, with accretion rates .

The magnitude of the planet mass that leads to gap opening and potential DZIB retreat is also of crucial importance for the IOPF theory. Gap opening has been studied by Lin & Papaloizou (1986, 1993), who proposed a “viscous-thermal” criterion that sets a gap opening mass. This mass is sensitive to the local disk viscosity. Magnetohydrodynamic simulations with realistic turbulence have also been carried out to study gap opening in turbulent disks (Nelson & Papaloizou, 2003; Zhu et al., 2013).

The second goal of this paper is to study gap opening and quantify the gap opening mass in the context of the IOPF disk model, i.e., an inner region of an active disk where the effective viscosity is undergoing a sharp increase in the radially inward direction. These gap opening masses can then be compared to masses of the innermost, so-called “Vulcan” planets. An initial comparison of the simple, analytic gap opening mass prediction with the observed Vulcans has been carried out by Chatterjee & Tan (2015, hereafter CT15, Paper II), who found a predicted scaling and normalization of planet mass versus orbital radius, , of , where is a dimensionless parameter that indicates the fraction of the Lin-Papaloizou mass scale that is assumed to lead to a deep enough gap to truncate planetary accretion of pebbles and initiate DZIB retreat. The normalization of this gap opening planet mass at given radius also depends on the Shakura-Sunyaev viscosity parameter, , in the dead zone inner boundary region, with the fiducial value of being adopted in the above formula. However, this value is quite uncertain and the observed planet masses may require a somewhat lower value of (Paper II), or a lower value of .

Finally, our third goal is investigate how gap opening may induce DZIB retreat, which then sets the location of second planet formation. Here we will present a simple, heuristic first exploration of this process by modeling the location of the DZIB as set by penetration of X-rays emitted from the protostellar corona to the disk midplane beyond the planet-induced gap. This location then sets the radius for an imposed radial profile of viscosity that simulates the transition region from an MRI-active inner region to an outer dead zone. Increased viscosity leads to reduced densities that make it easier for X-rays to penetrate further. We present example toy models of this process in which the location of the DZIB ends up stabilizing several initial gap widths away from the planet. Such separations can be compared to the observed separation of orbits of innermost STIPs planets.

In §2, we describe an analytic model to calculate our disk parameters. In §3, we describe our numerical set-up and test the numerical simulations against the analytic model. In §4, we study migration of planets of various masses located in the DZIB region. In §5, we study gap opening and DZIB retreat. In §6, we discuss the implications of our results for observed planets. We conclude in §7.

2. Analytic Estimates of Disk Structure and Gap Opening Mass

As in Paper I, we follow Frank et al. (2002)’s derivation of the structure of a viscously heated accretion disk. These results will be compared to the numerical simulations with FARGO, described below in §3. To achieve consistency with these simulations we find we need to make a small correction in the choice of the vertical optical depth equation and now adopt:


with the factor of 0.5 being introduced since the disk has two faces (or equivalently a change in the definition of “midplane” conditions). Thus when calculating the balance between energy dissipation and radiative cooling, the optical depth is integrated from the midplane to each surface of the disk. This change then leads to modest () changes in the normalization of the disk structure equations compared to CT14. For example, for disk mass surface density we now derive:


(CT14 derived a normalization of ), where is the mean particle mass (assuming ), is Boltzmann’s constant, is the power law exponent of the barotropic equation of state where we have normalized for with rotational modes excited, is Stefan-Boltzmann’s constant, is the stellar mass, is disk opacity (normalized to expected protoplanetary disk values, e.g., Wood et al., 2002), , (where is stellar radius), and is the accretion rate.

For the disk aspect ratio we find


The above equations will be used to set the initial conditions of the FARGO simulations described below.

Given this disk structure, the viscous criterion for the gap opening planet mass is :


Implementing our disk model, we obtain:


This is a factor smaller than the CT14 result.

The full set of disk structure equations with our revised normalizations are presented in Appendix A.

3. Numerical Simulation Methodology

We used the 2D code FARGO-ADSG (Baruteau & Masset, 2008a, b), which is built on FARGO (Masset, 2000), but with the energy equation and disk self-gravity implemented. Self-gravity has not been turned on in our calculation since the disk is far from gravitational instability in our problem, which focusses on inner regions near the star. We also implemented a simple method of radiative cooling that uses the disk surface density to calculate the optical depth and cooling rate as described in Zhu et al. (2012). We set up a disk with the same mass surface density profile, flaring index, opacity, adiabatic index, mean molecular mass and central stellar mass as the fiducial analytic model in §2.

The energy equation implemented in FARGO-ADSG is (Baruteau & Masset, 2008a):


where is the thermal energy per unit area, is the flow velocity and () denote heating (cooling) source terms, assumed to be positive quantities. In a steady accretion disk, these terms can be written as


where are components of the viscous stress tensor,

is the midplane temperature. Note is the optical depth from disk midplane to the surface7.

First we set up an accretion disk with a constant value of , i.e., without a transition zone, and a steady accretion rate of , simulating a range of radii from 0.02 to 0.3 AU. In all simulations we use the EVANESCENT boundary condition. This damps the disk values (surface density, velocities and energy density) to the initial axisymmetric disk conditions. The damping regions are rings within the radial ranges [] and [], where is the inner (outer) radius of the disk. The damping parameter is increased from zero to a maximum from the edge of damping zone to the edge of the disk.

Within each orbit, the damping amplitude is not a constant: the actual damping is the calculated damping function times a coefficient, which gradually grows from zero to one, linearly with time.

For our standard “medium resolution (MR),” the disk is evenly divided by 300 sectors azimuthally and logarithmically divided by 128 sectors radially. This gives a typical grid size of  AU, which is about the same as the Hill radius of a 1--mass planet located at 0.1 AU. We also carry out “high resolution (HR)” simulations that have twice the MR resolution, and a few test runs (of more limited time duration) at “super-high resolution” (SHR) that have four times the MR resolution.

We next set up a disk that has a radial jump in . Moving inwards, rises from at 0.1 AU to at 0.07 AU, i.e., over a transition width of  AU, with the transition described by part of the function from (at 0.1 AU) to (at 0.07 AU). The initial mass surface density profile is also chosen to give a steady accretion rate across this transition region. The steady radial profile achieved after 1000 orbits is close to this initial choice and is shown in Fig. 1.

Figure 1.— Mass surface density profile after 1000 orbits of our fiducial “transition” disk where rises from 0.001 at 0.1 AU to 0.01 at 0.07 AU with steady accretion rate of . The dashed line is the profile of a “constant ” disk (i.e., ).

4. Protoplanet Migration

We study the migration of protoplanets of various fixed masses (0.1, 0.5 and 1.0 ) that are inserted into the transition zone region of the disk by measuring the torques exerted on the planet by the disk. For each planet mass, we run simulations holding the planet at a constant radius, exploring uniformly from 0.085 to 0.115 AU with a spacing of 0.001 AU. So for each planet mass there are 31 MR and HR simulations run, for a total of 186 simulations. Note we keep the viscosity radial profile constant in all the simulations in this section. We run each simulation for 400 orbits, thus allowing the gas disk to achieve a quasi equilibrium structure. We then evaluate the torque on the planet by averaging over the next 100 orbits. Given the 2D nature of these simulations, we adopt a smoothing length in the calculation of the planet’s gravitational potential of (Müller et al, 2012).

The torque, , is scaled as , where is the orbital radius of the planet and is the average gas mass surface density at this orbital radius. The grid size in our standard MR run is more than twice the Hill radius of a planet with . We therefore carry out HR simulations with twice greater resolution and examine numerical convergence.

Figure 2.— Transition zone migration torque, , profiles for 0.1  (top), 0.5  (middle) and 1.0  (bottom) planets held at fixed radius. Medium resolution (MR) and high resolution (HR) runs are shown, with the torque induced by the disk on the planet being the average from orbit 400 to 500 after introduction of the planet. A positive torque leads to outward migration; a negative torque to inward migration. The HR results show decreasing from positive to negative values in the transition region (where starts decreasing as r increases just interior of 0.1 AU) . We expect the planet will migrate inwards from positions where the total torque is negative and migrate outwards where the total torque is positive. It will stop at zero torque. Thus these torque profiles are indicative of there being a stable “planet trap” in this location, where .

In all cases for 0.1 , 0.5  and 1.0 -mass planets (which correspond to very different levels of gap opening—see below), the torque profiles show a common behavior of decreasing from positive to negative values in the transition region (where starts decreasing as increases just interior of  AU), indicative of a stable “planet trap” in this location, where (see Figure 2). We note that the magnitude of decreases significantly as increases from the Type I to Type II regimes.

This confirmation of planet trapping at the DZIB transition zone is a key requirement of IOPF, since this allows a planet to continue to grow by pebble accretion from low to relatively high masses. Next we investigate how pebble accretion may be disrupted by gap opening by looking at the detailed structure of the gaps.

5. Gap Opening

5.1. Gap Structure & Opening Mass for Fixed Profile

The gap opening process in IOPF is critical to termination of pebble accretion of the first planet, thus setting its mass. It may also be important, in combination with dead zone retreat, in setting the location of the pressure maximum that leads to new pebble ring formation and then second planet formation.

We now investigate the disk structure that is induced by introducing planets of various masses at a fixed location of  AU, i.e., at the transition zone of the disk with a fixed profile. Figure 3 shows the radial profiles of mass surface density, midplane pressure and midplane temperature after 1000 orbits of evolution from introduction of planets of masses to in steps of 0.1 . We note that the gas profile settles very quickly, within orbits, to a profile close to that of the final state at 1000 orbits.

As planet mass increases, the radial profiles show a gradual deepening of the gap in both the mass surface density and pressure profiles. The temperature also dips due to reduced viscous heating and smaller optical depths. For this particular set-up, with a transition zone width of 0.03 AU, it is only once planet masses are that there is significant displacement of the azimuthally averaged pressure maximum away from the planet’s orbital radius. Table 1 lists the separation of the pressure maximum from the planet’s orbital radius at 0.1 AU, , but normalized by the planet’s Hill radius, . We define . As increases from 0.4 to 0.5 , grows by a factor of 10 from to . Figure 3 shows this also corresponds to the mass surface density maximum retreating outwards by a similar amount.

We have compared gap opening results for simulations with at MR and HR resolutions. We find very similar radial profiles of and , with agreement in values at better than 10% across the width of the gap region.

Figure 3.— Radial structure of disk mass surface density (top), midplane pressure (middle) and midplane temperature (bottom) during gap opening for planets of mass to in steps of 0.1  held fixed at the transition zone radius of 0.1 AU. Profiles are shown for medium resolution runs after 1000 orbits. Gap opening manifests itself via a decrease in , and in the vicinity of the planet and, once , a significant retreat of the local surface density and pressure maxima outwards by several Hill radii.

Figure 4 shows the 2D view of the perturbation of the disk mass surface density (compared to the disk with no planet) in the - plane. The strengthening spiral arms, along with deepening gap, are evident.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.61 0.48 0.42 0.38 4.75 4.47 5.09 4.86 4.68 5.27
Table 1Retreat of the local pressure maximum during gap opening
Figure 4.— 2D - view of the mass surface density perturbation, , induced by different planet masses from 0.1 to 1.0  (where is the surface density in the absence of a planet).

5.2. DZIB Retreat with an Evolving Profile

When a planet opens a gap in the disk leading to decreasing density in the vicinity of the planet, we expect higher ionization levels given the dependence of recombination rates, and thus greater likelihood of activating the MRI, which then results in the retreat of the dead zone inner boundary. Higher levels of ionization may also arise due to shock heating in spiral arms that are induced by the planet. We note that these spiral arms, propagating outwards into the dead zone, will by themselves also provide a source of enhanced effective viscosity. The lower densities of the gap region may also allow for increased penetration of X-rays that raise the ionization level and activate the MRI. Such effects could cause the viscosity to make its radial transition further out in the disk, which would lower the density in the gap region and vicinity even more.

The full treatment of these processes is a very complex problem, beyond the scope of the present paper. Here we present a simple, heuristic model for DZIB retreat based on X-ray penetration from the protostellar corona to the disk midplane. The X-rays are assumed to propagate from a height of 0.05 AU at (i.e., , which is expected to be –4 stellar radii for a solar mass pre-main sequence star). The thickness of the transition zone from the MRI-active region to dead zone is assumed to be characterized by , i.e., a mass surface density that scales with the X-ray penetration mass surface density due to absorption, .

Starting from the FARGO simulation with after 300 orbits (i.e., when a wide, but relatively shallow, gap has just been opened), the radius of the DZIB, , is re-evaluated each orbit by the condition that the mass surface density along the path from the X-ray corona to the midplane equals . For simplicity, we assume there is negligible material along this path length at locations interior to the planet, i.e., at  AU, which we consider to be a good approximation due to both viscous clearing of the inner disk and also since the disk’s vertical scale height is small in this region. Thus our path integral starts from  AU and extends out to . Note, that this path starts with a relatively shallow angle to the midplane of , which then decreases as the disk evolves and increases.

In addition, we need to specify the structure of the transition zone in the X-ray penetration region. We expect it to be broader than the transition zone set by thermal ionization of alkali metals. For simplicity, we assume the functional shape of the radial variation of is the same as that adopted in §3. The outer radius where starts rising from the DZIB value of 0.001 is set by the mass surface density along the pathlength being equal to . We note that the width of this transition region may be somewhat larger than the actual penetration depth of the X-rays (i.e., ), because of, for example, the outward propagation of turbulence, which has been seen in the simulations of Dzyurkevich et al. (2010). For the inner radius, , where reaches the full MRI value of 0.01, we investigate three cases, i.e., inner radii of 0.07 AU (Case A: the same as adopted for our migration and gap opening studies, above), 0.1 AU (Case B: the location of planet), and 0.11 AU (Case C: the location of the initial gap outer edge that has been induced by the planet).

For a simple estimate on the values of and to be used in this model, we note that X-ray opacity, even at  keV, may be dominated by heavy elements in the gas, rather than in dust, especially if dust has mostly coagulated into large grains and pebbles that have settled to the midplane and been trapped at the DZIB pressure maximum. For interstellar gas with solar abundances and the case of heavy elements depletion, Igea & Glassgold (1999) find an absorption opacity of (valid for –30 keV). Thus for requires and requires . Thus we choose a characteristic value of , equivalent to absorption optical depths of a few for  keV X-rays. We note, however, that the radiative transfer of these higher energy X-rays is likely to be controlled by Compton scattering by electrons present in neutral H, and He, which may allow increased penetration to the disk midplane at greater radial distances via scattering from the disk atmosphere. For this reason, along with the extra thickness of the transition zone that is expected due to propagation of turbulent wakes from the MRI-active zone (including vertically from modest scale heights), we choose a fiducial value of .

From the above discussion, it is clear that more detailed calculations are needed to better constrain the combined parameter , including ionization equilibrium calculations, e.g., similar to those of Igea & Glassgold (1999) and then applying the resulting resistivities to determine the extent of the MRI-active zone (Mohanty & Tan, in prep.). The time variable nature of the X-ray luminosity may also need to be accounted for.

Figure 5.— Radius of DZIB (i.e., the outer transition radius of profile) as set by X-ray penetration in evolving disks. Medium resolution (MR) results are shown with dashed lines; high resolution (HR) results with solid lines. The three cases of inner transition radii of 0.07 AU (Case A, red lines), 0.1 AU (Case B, blue lines) and 0.11 AU (Case C, green lines) are shown. In each case, the estimated converged, final retreat radius is shown with a horizontal dot-dashed line (see text).

Figure 5 shows results from these DZIB retreat models. We find that increases gradually over timescales of several thousand orbits of the planet. For Case A, by 6500 orbits the rate of increase is very slow, i.e., from the 6500th to 7000th orbit the percentage increase is 0.06% for the MR run and 0.07% for the HR run. The equivalent results for Cases B and C are similar. However, the value of after 7000 orbits does depend somewhat on the resolution of the simulation: for the Case A MR run  AU; for the HR run  AU, i.e., the amount of retreat () is about 38.6% greater in the HR run. We thus also carried out twice higher, SHR, resolution runs for 1000 orbits. We find that is 5% greater than in the HR simulation. We use these results to estimate that the final numerically converged retreat distance is 5-10% greater than the HR result (or about 40% greater than the MR result), i.e.,  AU, i.e.,  AU. For this case with , for which the planet’s Hill radius is  AU then , which is a significantly greater retreat (i.e., greater) than due to simple gap opening in the fixed model.

For Cases B and C the amount of retreat is larger and takes somewhat longer to stabilize (although the retreat has essentially stopped by 7000 orbits). After 7000 orbits, the HR simulations find a 24.0% and 20.8% larger retreat than the MR runs for Cases B and C, respectively. The SHR results at orbits are 5.6% and 4.9% greater, respectively. Thus we estimate final converged values are about 25% greater than the MR results, i.e., for Case B DZIB retreat to 0.20 AU ( AU; ) and for Case C DZIB retreat to 0.22 AU ( AU; ). These results are also summarized in Table 2.

As a check on self-consistency, we have also investigated the potential effect of X-ray penetration inducing an evolving profile in a disk that does not yet have a planet or gap. For this model, where the accretion disk extends uninterrupted to very close to the protostar at  AU, we set the height of the stellar X-ray corona to be equal to this inner disk radius, i.e., 0.02 AU. We also start the path integral for X-ray penetration from this distance. We find that for the same choices of and adopted above, the X-rays do not penetrate beyond  AU and so the DZIB profile set by thermal ionization of alkali metals is not expected to be affected. However, we find that utilizing a larger X-ray corona height of 0.05 AU does enable X-ray penetration to beyond 0.1 AU in the no-planet disk. Thus, in the context of this simple model, DZIB retreat is also influenced by properties of the scaleheight of X-ray emission from the protostar.

Overall, we therefore see that the amount of DZIB retreat, which in the IOPF scenario is expected to set the location of second planet formation, depends somewhat sensitively on the detailed structure of the transition zone and the X-ray emission properties of the protostar. Thus future work to model this zone including the full physics of activation of the MRI in response to a changing ionization fraction, e.g., due to X-rays or other ionization sources, is needed. However, we can still conclude that, for the parameters explored here, the amount of retreat is significantly greater than the retreat due to initial gap opening. In §6 we compare these retreat distances to the orbital spacings of the Kepler-observed planets.

Medium Res.(7000 orbits) High Res.(7000 orbits) Estimated Convergence
Case A (, =0.07 AU) 15.9 (0.028 AU) 22.1 (0.039 AU) 25.3 (0.045 AU)
Case B (, =0.10 AU) 41.2 (0.073 AU) 51.1 (0.091 AU) 56.3 (0.100 AU)
Case C (, =0.11 AU) 52.4 (0.093 AU) 63.3 (0.112 AU) 67.6 (0.120 AU)
Table 2DZIB retreat for evolving models for different inner transition radii,

Figure 6 shows the radial structure of the disks for these evolving models after 7000 orbits, and compares to the fixed profile model, which was already settled after 1000 orbits. This figure shows that the local pressure maximum remains closely located with the outer transition radius, i.e., where it begins to rise from its DZIB value of 0.001.

Figure 6.— Radial structure (Left panel: mass surface density; Right panel: midplane pressure) of a disk with a 0.5 planet held fixed at 0.1 AU. 1st row: The fixed profile disk model ( is shown by the dotted red line). Disk structure is shown after 1000 orbits (solid blue line), compared to initial state (dashed green line). The extent of the planet’s Hill sphere is shown by the vertical shaded band. 2nd row: The evolving disk model (HR run) at 7000 orbits due to X-ray penetration in which the inner transition radius is held fixed at 0.07 AU (Case A). 3rd row: Same as above, but for Case B with inner transition radius at 0.1 AU. Bottom row: Same as above, but for Case C with inner transition radius at 0.11 AU.

Figure 7 shows the 2D - plots of perturbations of mass surface density and pressure for the Case A evolving model and compares to the fixed model, both of which have . Figure 8 shows the absolute values of these quantities. For completeness, to better illustrate physical appearance of disk structures, we also show true spatial maps of disk mass surface density structure of Case A in Figure 9. These figures reveal the decreasing absolute mass surface density of the induced spiral arm in the case where DZIB retreat has proceeded further.

Figure 7.— Mass surface density (top row) and midplane pressure (bottom row) perturbation plots of disks with a planet with (Left column) fixed profile (after 1000 orbits) and (Right column) evolving profile (Case A, after 7000 orbits).
Figure 8.— Same as Figure 7, but now showing absolute values of mass surface density (top row) and midplane pressure (bottom row) of disks with a planet with (Left column) fixed profile (after 1000 orbits) and (Right column) evolving profile (Case A, after 7000 orbits).

We emphasize that we have presented a very simple, heuristic model for a disk in which the viscosity profile responds actively to X-ray penetration (which itself is sensitive to the disk’s structure). More accurate modeling would involve self-consistent calculation of the viscosity caused by the MRI in the transition zone (e.g., Bai & Stone, 2011; Mohanty et al., 2013), including calculation of the ionization structure of the disk as it is affected by gap opening. The results that are needed from such modeling are the shape and absolute values of the disk viscosity in the MRI transition zone, which could then be implemented parametrically in, e.g., FARGO hydrodynamic simulations of the disk. However, potential instability of the dead zone inner edge, which may lead it to vary its radial location (Latter & Balbus, 2012), or vortex formation at the inner edge (Faure et al., 2015), should also be assessed. Although the model presented so far should be regarded as illustrative, as we discuss in §6, the scale of pressure maximum retreat that it exhibits is comparable to the orbital separations of the innermost planets in STIPs.

Figure 9.— Same as Figure 8 top row, but now showing a spatial map of the disk structure induced by a planet with (Left) fixed profile (after 1000 orbits) and (Right) evolving profile (Case A, after 7000 orbits).

6. Implications for Observed Planets

CT14 examined the separation of planet orbits in STIPs (normalized by Hill radius of inner planet in the pair), finding a broad range of separations from innermost, Vulcan, planet to next planet of few–100, with a broad peak from . Note, however, that this analysis was based on the simple, single power law planetary mass-size relation of Lissauer et al. (2011).

CT15 analyzed the available mass-size data for STIPs planets and derived an improved, piece-wise power law mass-size relation. Here, we use this relation (PL3 in CT15) to re-evaluate the distributions of through , where the numerical value of the subscript refers to the number of the planet that is the innermost in the pair being considered, counting out from the star. These results are presented in Figure 10.

Note some of the dispersion in separations is likely to be induced by an intrinsic dispersion in the densities of the planets, leading to inaccurate estimation of masses. For example, the observed dispersion in density of a factor of about 5 (e.g., CT15) leads to a dispersion in inferred mass of the same factor and a dispersion in and of a factor of 1.7. Also, the observed values of may be overestimated if some fraction of second planets are non-transiting and thus not detected due to slight orbital misalignments.

In Figure 10 we also show the locations of for the fixed and evolving disk models we considered. We note that the evolving disk models can achieve separations of to 70 , which overlap with the observed distribution of .

As discussed by CT14, the distributions of are significantly different from those of the outer separations. This conclusion remains unchanged by our use of the improved piecewise power-law CT15 mass-size relation. For KPC systems with , the probability that and are drawn from the same distribution is . Similarly, for systems with , the respective probabilities that two distributions are drawn from the same underlying distribution are , , and for and , and , and and . This indicates that the distribution of separations between the first two planets is different than between the 2nd and 3rd and other outer pairs of planets in these systems.

As also discussed by CT14, the generally larger values of may be a signature of the greater relative effect of inner disk clearing after first, Vulcan planet formation on the location of the second planet, compared to later, more incremental dead zone retreats.

Figure 10.— A comparison of orbital separations of the Kepler-detected planets (see legend) and our simulation results for DZIB pressure maximum retreat: vertical long-dashed line shows for the initial gap-openning in the fixed profile disk. The three short dashed lines show the estimated values of for Cases A, B, C, respectively, of the evolving profile disk.

7. Summary and Discussion

We have explored three main aspects of the Inside-Out Planet Formation (IOPF) scenario with numerical simulations of planet-disk interactions. Our main findings are the following:

(1) Planets of all relevant masses are trapped at the dead zone inner boundary transition region, so planetary migration (of Type I or Type II) does not appear to be a problem for IOPF. Especially during the early stages of potential Type I migration, the fact that the protoplanet is trapped at its formation location should allow it to continue to accrete, especially by pebble accretion.

(2) The mass scale of significant gap opening that can affect pebble accretion by displacing the local pressure maximum away from the planet is , i.e., , which is about in the fiducial case at . This is similar to the fiducial value adopted by CT14, especially since given the renormalization of by a factor of 0.745 in this paper (§2), this means that the overall change in planet mass is only a factor of 1.24. However, these results are derived in the context of pure hydrodynamic simulations and may be sensitive to including extra physics of MRI activation in the gap region.

(3) A simple model for MRI activation in the gap region, implemented by an evolving viscosity profile set by X-ray penetration, leads to greater retreat of the pressure maximum out to separations of to 70 Hill radii of the planet, but sensitive to model assumptions about the shape of the viscosity profile in the MRI-active region to dead zone transition region. Such separations overlap with the observed orbital separations of innermost planets in STIPs. Future work is needed to better model viscosity profiles in such transition regions.

While the above results are supportive of the IOPF model for being relevant for formation of STIPs, a number of open questions remain. For example, even though planets with masses appear likely to be trapped at a location where they can continue to grow by pebble accretion, it is possible that this is not the case at earlier stages for lower mass protoplanets. Faure et al. (2015) have pointed out that inwardly migrating vortices can form at the dead zone inner edge and they can potentially interact with the planet, possibly causing it to also migrate inwards. However, their 3-D simulation results are for an unstratified disk and the mass accretion rate is not a constant across the dead zone inner edge, leading to mass pile up and potentially influencing vortex formation. The potential effects of vortex interactions still needs to be confirmed in 3-D stratified simulations where the presence of the surface active layer may transport the disk mass smoothly across the dead zone inner edge. In our 2-D simulations, we used a small viscosity in the dead zone region to represent the effects of an active layer and/or turbulent wakes spreading from the inner MRI-active region, which may control the accretion rate in this radial region of the disk. With this set-up we do not observe the sharp density peaks at the dead zone inner edge that were present in the Faure et al. (2015) simulations.

Further work is needed on the IOPF model to study the transition of the pebble ring to a single dominant protoplanet, potentially involving streaming and gravitational instabilities of the pebble population and/or Rossby wave instabilities in the gas, followed by oligarchic growth of planetesimals combined with continued pebble accretion. Improved study of pebble supply truncation during the process of gap opening, including utilizing 3-D simulations (see Lambrechts et al., 2014, for example calculations focussed on the outer disk), are also needed.

Potential migration of planets after their main accretion phase and retreat of their natal DZIB, including due to interactions between neighboring planets, is another topic for future study. However, we note that for the fiducial IOPF model, the gas mass that is present in the vicinity of first planet to form is relatively small compared to the planet’s mass, thus limiting the scope of its migration.

We thank F. Masset and S. Mohanty for helpful discussions, as well as the comments of an anonymous referee. We acknowledge support from NASA ATP grant NNX15AK20G (PI: JCT). The authors acknowledge University of Florida Research Computing for providing computational resources and support that have contributed to the research results reported in this publication.

Appendix A Derivation of disk properties

Here we update the analysis presented in CT14, which involves modest changes in disk properties, as discussed in §2. Conservation of angular momentum in a viscous disk implies:


where is the shear torque from viscosity:


In a steady disk, we have , so integrating A1 we obtain:


Substituting , we obtain:


The boundary condition at , , so constant is:




For , we have , and substituting into A4, we have:


So the energy dissipation rate per unit area (only via one face of the disk) at radius r is:


From the mid-plane to disk surface, the optical depth is:


From Eq. 3.11 in Hubeny (1990), we know the relation between surface temperature and actual temperature where the optical depth is to the surface of the disk:


Here , and in our situation, , , so


We have an expression for mid-plane temperature:


but we still have an unknown profile, so now we try to deduce it. In an disk, we have:


Substituting in A7, we obtain:


For an adiabatic sound speed , we have:


Then the mid-plane pressure satisfies:


Now we have the following relation between , and :


so, substituting for from equation A12, the surface density profile is:


and mid-plane temperature (here we use instead of ) is:


For scale height , usually it is assumed that the disk is vertically isothermal in a passive disk, here in an active disk, we use an adiabatic vertical structure:


so the aspect ratio is:


Since the disk vertical structure satisfies , mid-plane density can be expressed as:



  1. affiliation: Department of Astronomy, University of Florida, Gainesville, FL 32611
  2. affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544
  3. affiliation: Department of Astronomy, University of Florida, Gainesville, FL 32611
  4. affiliation: Department of Physics, University of Florida, Gainesville, FL 32611
  5. affiliation: Department of Astronomy, University of Florida, Gainesville, FL 32611
  6. affiliation: Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Physics and Astronomy, Northwestern University, Evanston, IL 60208
  7. During numerical testing we found an error in the public version of the code: it only calculates half the value of .
  8. footnotetext: The separation between new pressure maximum and planet orbital radius, scaled to planet’s Hill radius. These results are for disks with fixed profile with a transition width of  AU. Note, location of pressure maximum is measured at the center of grid cell of maximum pressure.
  9. footnotetext: measures the separation between the planet (at 0.1 AU) and the outer transition of profile (as set by X-ray penetration) in units of the planet’s Hill radius. The outer transition radius is at a location very similar to that of the local pressure maximum (see Fig. 6).


  1. Bai, X.-N. & Stone, J. M.2011, ApJ, 736, 144B
  2. Baruteau, C. & Masset, F. 2008, ApJ, 672, 1054B
  3. Baruteau, C. & Masset, F. 2008, ApJ, 678, 483B
  4. Baruteau, C., Crida, A., Paardekooper, S. J. et al. 2014, Protostars and Planets VI, Univ. Arizona Press, Tucson, p.667
  5. Bitsch, B., Morbidelli, A., Lega, E., Kretke, K. & Crida, A. 2014, A&A, 570A, 75B
  6. Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582A.112B
  7. Chatterjee, S., & Ford, E.  B. 2015, ApJ, 803, 33C
  8. Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53
  9. Chatterjee, S., & Tan, J. C. 2015, ApJ, 798, L32
  10. Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444
  11. Cossou, C., Raymond, S. N. & Pierens, A. 2013, A&A, 553L, 2C
  12. Cossou, C., Raymond, S. N. & Pierens, A. 2014, IAUS, 299, 360C
  13. Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70
  14. Fabrycky, D. C., Lissauer, J. J., Ragozzine, D. et al. 2014, ApJ, 790, 146
  15. Fang, J. & Margot, J. L. 2012, ApJ, 761, 92
  16. Faure, J., Fromang, S., Latter, H. & Meheut, H. 2015, A&A, 573A, 132F
  17. Frank, J., King, A. & Raine, D. J., 2002, Accretion Power in Astrophysics: Third Edition (Cambridge University Press; ISBN 0521620538)
  18. Goldreich, P. & Schlichting, H. E. 2014, AJ, 147, 32G
  19. Hansen, B. & Murray, N. 2013, ApJ, 775, 53
  20. Hansen, B. & Murray, N. 2012, ApJ, 751, 158
  21. Hubeny, I. 1990, ApJ, 351, 632H
  22. Igea, J., & Glassgold, A. E.1999, ApJ, 518, 848
  23. Kley, W. & Nelson, R. P. 2012, ARA&A, 50, 211
  24. Kretke, K. A. & Lin, D. N. C. 2012, ApJ, 755, 74K
  25. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107
  26. Lambrechts, M. & Johansen, A. and Morbidelli, A. 2014, A&A, 572A, 35L
  27. Latter, H. N. & Balbus, S. 2012, MNRAS, 424, 1977L
  28. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322
  29. Lin, D. N. C., & Papaloizou, J. C. B. 1986, ApJ, 307, 395
  30. Lin, D. N. C., & Papaloizou, J. C. B. 1993, Protostars and Planets III, Univ. Arizona Press, Tucson, p. 749
  31. Lissauer, J. J., Ragozzine, D., Fabrycky, D. C. et al. 2011, ApJS, 197, 8L
  32. Masset, F. S. 2000, A&A, 141,165
  33. Masset, F. S. & Morbidelli, A. and Crida, A. and Ferreira, J. 2006, ApJ, 642, 478
  34. Matsumura, S., Pudritz, R. E. & Thommes, E. W. 2009, ApJ, 691, 1764
  35. Mohanty, S., Ercolano, B., & Turner, N. J. 2013, ApJ, 764, 65
  36. Mullally, F., Coughlin, J. L., Thompson, S. E. et al. 2015, ApJS, 217, 31M
  37. Müller, T. W. A., Kley, W. & Meru, F. 2012, A&A, 541,123
  38. Nelson, R. P. & Papaloizou,J. C. B. 2003, MNRAS, 339, 993N
  39. Ogihara, M., Morbidelli, A. & Guillot, T. 2015, A&A, 578A, 36O
  40. Raymond, S. N. & Cossou, C. 2014, MNRAS, 440L, 11R
  41. Schlichting, H. E. 2014, ApJ, 795L, 15S
  42. Toomre, A. 1964, ApJ, 139, 1217T
  43. Varnière, P. and Tagger, M. 2006, A&A, 446L, 13V
  44. Wood, K., Wolff, M. J., Bjorkman, J. E., & Whitney, B. 2002, ApJ, 564, 887
  45. Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459Y
  46. Zhang, X., Liu, B., Lin, D. N. C. & Li, H. 2014, ApJ, 797, 20Z,
  47. Zhu, Z., Hartmann, L., Nelson, R. P. & Gammie, C. F. 2012, ApJ, 746, 110Z
  48. Zhu, Z., Stone, J. M., & Rafikov, R. R. 2013, ApJ, 768, 143
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
The feedback must be of minumum 40 characters
Add comment
Loading ...

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description