Pebble Evolution and Formation Timescales

Inside-Out Planet Formation. IV. Pebble Evolution and Planet Formation Timescales


Systems with tightly-packed inner planets (STIPs) are very common. Chatterjee & Tan (2014) proposed Inside-Out Planet Formation (IOPF), an in situ formation theory, to explain these planets. IOPF involves sequential planet formation from pebble-rich rings that are fed from the outer disk and trapped at the pressure maximum associated with the dead zone inner boundary (DZIB). Planet masses are set by their ability to open a gap and cause the DZIB to retreat outwards. We present models for the disk density and temperature structures that are relevant to the conditions of IOPF. For a wide range of DZIB conditions, we evaluate the gap opening masses of planets in these disks that are expected to lead to truncation of pebble accretion onto the forming planet. We then consider the evolution of dust and pebbles in the disk, estimating that pebbles typically grow to sizes of a few cm during their radial drift from several tens of AU to the inner, AU-scale disk. A large fraction of the accretion flux of solids is expected to be in such pebbles. This allows us to estimate the timescales for individual planet formation and entire planetary system formation in the IOPF scenario. We find that to produce realistic STIPs within reasonable timescales similar to disk lifetimes requires disk accretion rates of and relatively low viscosity conditions in the DZIB region, i.e., Shakura-Sunyaev parameter of .

Subject headings:
protoplanetary disks, planets and satellites: formation


1. Introduction

More than 4000 exoplanet candidates (e.g., Mullally et al., 2015) have been discovered by Kepler since 2009, a significant percentage () of which are in systems with tightly-packed inner planets (hereafter STIPs). STIPs are very different from our Solar System: they typically contain three or more detected planets of radii and with periods less than 100 days (Fang & Margot, 2012), i.e., with orbital radii of  AU.

One way to form STIPs may be via inward migration of planets formed in the outer disk (e.g., Kley & Nelson, 2012; Cossou et al., 2013, 2014). However, it can be quite difficult to concentrate planets close to the host star to the degree of observed STIPs (McNeil & Nelson, 2010). Also migrating planets tend to become 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). This lack of resonant pile-up has motivated studies of mechanisms that cause a lower efficiency of resonance trapping (Goldreich & Schlichting, 2014) or breaking of resonance by later dynamical interaction with, e.g., a planetesimal disk (Chatterjee & Ford, 2015), by the unstable nature of compact resonant chains (Izidoro et al., 2017), or by magnetospheric rebound from the host star during gas disk dispersal (Liu et al., 2017).

As an alternative to inward migration, in situ formation scenarios (Hansen & Murray, 2012, 2013) have considered a group of protoplanets that are initially distributed inside 1 AU with high concentration. After following dynamical evolution for 10 Myr, collisional oligarchic growth can lead to systems quite similar to STIPs. However, how this initial condition was established was not well studied in these works and may be subject to some difficulties, such as early triggering of gravitational instability (Raymond & Cossou, 2014; Schlichting, 2014) (see also, Chiang & Laughlin, 2013). In addition, simulations that account for the effects of gas on protoplanet migration during the oligarchic growth phase find that systems are produced with planet masses decreasing steeply with orbital radius, which is different from the relatively flat scaling seen in STIPs (Ogihara et al., 2015).

The Inside-Out Planet Formation (IOPF) scenario proposed by (Chatterjee & Tan, 2014, hereafter Paper I) is a new type of in situ formation model. It starts with pebble delivery to the disk midplane transition region between an innermost magnetorotational instability (MRI) active zone and a nonactive “dead zone,” where there is a local pressure maximum. The pebbles are trapped here and build up in a ring, which then forms a first 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, which may only be a partial gap with a modest reduction in gas mass surface density, pushes the pressure maximum outwards by at least a few Hill radii from the first planet, thus creating a new pebble trap, and thus location of second planet formation (Hu et al., 2016, hereafter Paper III) (see also Lambrechts et al., 2014). Furthermore, MRI-activation in this gap region beyond the first planet may induce enhanced viscosity and thus further outward retreat of the dead zone inner boundary (DZIB) pressure trap, leading to larger separations to the orbit of the second planet.

The pebble delivery rate is crucial for setting the timescale of pebble ring and planet formation in the IOPF scenario. To calculate this, we need to model the evolution of dust and pebbles in the protoplanetary disk and determine if the majority of solid material, which in the outer disk is initially in small dust grains, become incorporated into pebbles by the time it reaches the inner region. The IOPF model requires dust grains to grow into pebbles as smaller grains are well coupled to the gas and will not be trapped at a local pressure maximum (e.g., Zhu et al., 2012). The size distribution of pebbles is also important for models of planetesimal formation, e.g., via the streaming instability. This may be triggered in both the pebble ring as a first step for in situ planet formation, or in the outer disk as a potential first step for outer planet formation that may eventually truncate pebble delivery to the inner disk. Finally, modeling of dust and pebble evolution is a necessary step for then providing observational predictions of the emission from the disk during various phases of IOPF.

There have been a number of previous studies of dust and pebble evolution in global disk models, including radial drift and particle interaction (e.g., Brauer et al., 2008; Birnstiel et al., 2010). Birnstiel et al. (2012) proposed a two population model that divides dust into one smaller fixed size group and one larger variable size group that represents grain growth. Sato et al. (2015) proposed a simplified pebble-pebble interaction scenario by reducing the full Smoluchowski equation to a single sized evolution equation, yielding pebble fluxes and typical sizes in disks on scales down to 1 AU. A Lagrangian dust evolution model has been investigated by Krijt et al. (2016), which tracks batches of particles in the disk as they drift radially and grow through collisions.

In this paper, we present a new pebble-dust evolution model and apply it specifically to calculate the timescales needed to build up a STIP forming via IOPF, i.e., limited by the rate of pebble delivery. In §2, we review the equations that govern the disk models that are relevant to IOPF, including presenting new disk structure models that include realistic opacities and the transition from an “active,” accretion-powered inner disk to a “passive,” stellar-heated outer disk. In §3, we discuss a simple estimate of the timescale of first, so-called “Vulcan” planet formation in steady accretion rate disks, including a reanalysis of the gap-opening mass criterion that was studied in Paper III. We also estimate the size of the region of the disk that is needed to supply the material for this planet. In §4, we model the radial drift and sweep up growth of a single pebble in our model disks. In §5, we introduce our full pebble evolution model and numerical results for different disk parameters. In §6 we outline the assembly of a STIPs-like system with pebble flux and planet mass criterion derived from earlier sections. We discuss the implications of our results and summarize our conclusions in §7.

2. Disk Structure

The structure of a steady accretion disk can be specified with a given accretion rate , Shakura & Sunyaev (1973) viscosity parameter , stellar mass , and midplane temperature versus radius relation (in particular specifying whether the energy source of the disk is heating mostly by its own accretion, i.e., an “active” disk, or whether it is heating mostly by stellar irradiation, i.e., a “passive” disk):


where is the midplane sound speed, is the power law exponent of the barotropic equation of state , is Boltzmann’s constant, is the mean particle mass (with fiducial value set by assuming ), is the disk vertical scaleheight, is the local Keplerian speed, is the viscosity, is the Shakura-Sunyaev dimensionless viscosity parameter (with fiducial value set for dead zone conditions, but we will consider other values also), is the gas mass surface density, is the accretion rate (with fiducial normalization set by observed accretion rates of transition disks, but again we will consider a range of values; see Paper I) and is midplane gas density. Only equations (1) and (2) depend on the energy source and equation of state, while the rest are general equations followed by all viscous disks.

The solution for the structure of an active accretion disk is achieved by balancing viscous thermal dissipation with vertical radiative cooling (per unit area from one face of the disk):


This generally applies in the inner disk. In this case the temperature equation is (Paper I):


and gas mass surface density is given by:


where with fiducial normalization to a value of 1.4 for with rotational modes excited, is the disk mean opacity (with fiducial normalization here appropriate for inner disk conditions; but, note we will more generally use tabulated opacities from Zhu et al., 2009), is the Stefan-Boltzmann constant, is the stellar mass with fiducial normalization of , , is the stellar radius, and . The two parameters with some of the largest uncertainties are and , so we shall consider the effects of varying their fiducial values.

In the outer disk, stellar irradiation is expected to become the dominant energy source compared to local viscous accretion heating. Following Chiang & Goldreich (1997), the temperature structure of a stellar irradiated, passive disk follows , which is much shallower than the relation of the active disk.

To derive in the passive disk regime, the following equations are solved iteratively:


Here is the “grazing angle,” i.e., the angle of incidence of stellar radiation at the disk surface, is the surface, photospheric temperature of the star, is the flaring index, i.e., , and is the height of the disk photosphere above the midplane. Usually, is a good approximation in most parts of the disk (Chiang & Goldreich, 1997). If we choose a solar mass T-Tauri star with temperature and radius , we have the disk temperature versus radius equation


which is independent of and .

For our full numerical solutions, we set up a hybrid disk that is solved first in its inner region as an active disk, and then calculate a transition radius where the contribution from stellar irradiation becomes equal to that due to accretion heating. Beyond this transition radius, we scale the temperature in the manner appropriate for the passive disk regime, i.e., in particular, we compare the midplane temperatures of the active and passive disk models and choose the hotter one as the actual value.

The results of these constant (the fiducial value) and and constant , and disk models are shown in Figure 1. Such values of have been inferred in the inner regions of simulated dead zones (Dzyurkevich et al., 2010). Whether such dead zone conditions and values of continue to apply in larger scale regions of the disk remains an open question (e.g., Mohanty et al., 2013), however, for simplicity we will assume that they do and calculate the structure of the disk out to 30 AU. The range of accretion rates that we consider is guided by observations of transition disks (e.g., Manara et al., 2014). Another risk of applying constant and stable accreting disk to larger radius is overestimating the disk mass. In our fiducial case, as seen in Figure 1, the dust mass within 30 AU is about 100 , while recent observations on dust continuum (van der Marel et al., 2016) give an average around 40 for most transitional disks. This factor of 2 should be considered when analyzing the simulation result. The active to passive transition radius for the fiducial disk () is 1.6 AU. We note that the model assumption of disk vertical optical depth being breaks down in the cases with and in the outer disk, e.g., for this occurs for  AU. Note also that an “inner transition” in disk structure at 0.14 AU results from the drop in opacity due to dust sublimation. However, in the context of IOPF, we expect to also rise dramatically when K (this location is referred as hereafter) as the MRI is activated by thermal ionization of alkali metals. This change in has not been included in the models shown in Fig. 1 since our main focus here is modeling the region that is the supply reservoir for solid material delivered to 0.1 AU locations. Still, we have extended the disk models within 0.1 AU, since for the low accretion rate case (), .

In fact, there are some close-in planets with orbits within 0.1 AU, e.g., Kepler-10b (Dumusque et al., 2014) is a super-Earth planet orbiting a solar type star at 0.016 AU. For such planets to form in situ places constraints on IOPF models. For example, such systems would need to form in a disk with relatively low accretion rate that has a small value of (see below).

Figure 1.— Structural profiles of IOPF disk models used in this paper, assuming (left column) and (right column), for accretion rates (red dotted), (black solid) and (blue dashed). All models are for a one solar mass central star. From top to bottom, the rows show: gas mass surface density (); midplane temperature (); midplane opacity () (assumed to be vertically constant); disk aspect ratio (); enclosed mass in solids (), i.e., initially assumed to be dust, summed in the disk from out to radius for a solid to gas mass ratio of . Note that the model assumption of disk vertical optical depth being breaks down in the cases with and in the outer disk, e.g., for this occurs at  AU.

3. “Vulcan” planet formation

Our goal in this section is to obtain a simple estimate for the formation time and the size of the supply reservoir of the first, i.e., innermost planet to form in the IOPF scenario (Chatterjee & Tan, 2015, hereafter Paper II). To do this requires knowing the mass of such planets and the mass flux of pebbles. Here we will assume a steady accretion rate with most mass flux of solids being in the form of pebbles. The validity of this assumption will be explored in the following sections of the paper.

3.1. Location

The radial location of Vulcan planet formation is assumed to be set by where the temperature reaches 1,200 K that leads to thermal ionization of alkali metals, especially Na and K. This location varies depending on the accretion rate and other disk properties. In an active accretion disk, which is the relevant regime for the optically thick inner disk, this location is (Paper I):


where is a dimensionless parameter of order unity to account for potential differences from the pure viscous disk model, e.g., due to the effects of extra energy extraction on midplane temperature by a disk wind. Paper I and II adopted a fiducial value of of 0.5. In Paper III, we used 0.1 AU as the typical location of Vulcan planets for hydrodynamic simulations of gap opening, where , , which corresponds to .

3.2. Mass

The mass of the Vulcan planet (and all planets forming via IOPF) is assumed to be set by gap opening, which then leads to displacement of the local pressure maximum away from the planet, to a larger radius in the disk, and thus the truncation of pebble accretion. Paper III investigated this process with 2D hydrodynamic simulations for the case of the fiducial accretion rate (), a value of for the inner dead zone that then rises to in the MRI-active region, and for the planet set at a fixed location of 0.1 AU from a star. The mass of the planet leading to the first significant displacement of the pressure maximum was assessed relative to the viscous thermal criterion gap opening mass of Lin & Papaloizou (1993):


Here . The last two numerical evaluations of this mass illustrate its relatively strong dependence on . Paper III found that, for the case, the critical pressure-displacement gap opening mass was about 50% of the viscous-thermal gap opening mass (i.e., this fraction was expressed as ). Note that the mass evaluated in equation (14) is assessed only from the viscous criterion, but, for the conditions with , is consistent with the thermal criterion also, i.e., the Hill radius of the planet is similar in size to the disk scale height. However, as we see now, this condition breaks down when , as the planet mass is much lower.

To better capture the behavior of gap opening for more general conditions in which the thermal criterion also begins to limit gap opening, we now improve our treatment by adopting the expression of Duffell (2015) for the gap opening mass (same scalings of and is also presented in Zhu et al. (2013) and Fung et al. (2014)):


where is a dimensionless parameter of order unity (note this absorbs the dimensionless factor, , introduced by Duffell (2015)). Note that the scaling of this relation is slightly different from the pebble isolation mass in Lambrechts et al. (2014), which depends on , but not . However, a more recent prescription for the pebble isolation mass in Bitsch et al. (2018) also includes a dependence on . A major difference of our results from these other studies is that our planet gap opening mass measurements are done specifically at a region in the disk where there is a transition in the value of that causes a strong change in . On the other hand, the results of Lambrechts et al. (2014) and Bitsch et al. (2018) are for disks with constant and thus a relatively slowly varying profile.

To test equation (15) and constrain , we now extend the simulations of Paper III to investigate the gap opening mass for different values of and . We perform a series of FARGO (Masset, 2000) hydrodynamic simulations on disks with an viscosity transition. The inner region is still assumed to be an MRI-active zone with . Beyond this is the dead zone with a much lower value of . We implement the transition following the same methods as Paper III, with the transition set to be at 0.1 AU, which is the location of the local pressure maximum before the introduction of a planet in the disk. Then planets of different masses are inserted in the disk and held fixed at this location (note Paper III showed that growing planets would be trapped here given the torques exerted by the disk). We study the cylindrically averaged pressure profiles and identify the critical gap-opening mass for IOPF as corresponding to the point when the local pressure maximum in the disk becomes significantly displaced from the planet, i.e., by more than its Hill radius,


where is the planet mass normalized to one Earth mass.

Figure 2.— Steady state, cylindrically averaged disk midplane pressure profiles, starting with an initial case with no planet. The pressure maximum, due to the dead zone inner boundary, has been set at 0.1 AU. We set in the inner MRI-active disk, transitioning to in the dead zone. Then the effects of adding planets of increasing mass and held fixed at 0.1 AU are studied. These planets open gaps of increasing depth, which eventually leads to displacement of the local pressure maximum away from the planet. This occurs for , i.e., .

Figure 2 shows examples of these pressure profiles for the case of a disk with and a transition from to in the vicinity of 0.1 AU. We see in this case that a planet with mass of is able to open a deep enough gap that the local pressure maximum retreats significantly, i.e., by AU (in this case, Hill radii) from the planet. This mass scale corresponds to a value of in the normalization of the Duffell (2015) gap opening mass (eq. 15). In the IOPF model this means that the next planet may form from a new pebble ring at this location or somewhat further out if the MRI-active region also spreads outwards in the lower densities enabled by gap opening.

Figure 3a shows the results of four different series of simulations, including the example shown in Fig. 2, in which in the dead zone is varied from to and the gap opening masses (by the pressure maximum displacement criterion) are assessed and compared to those predicted by the simple viscous criterion (eq. 14) and the Duffell (2015) estimate (eq. 15). Figure 3b shows the results of a similar series of gap opening simulations, but now varying from to . We find the simple viscous criterion underestimates the planet mass required for displacing the pressure maximum in both the low and high regimes. These are conditions when the disk is relatively denser and hotter and its vertical scale height is relatively large compared to the Hill sphere of the planet. Our results are consistent with those of Duffell & MacFadyen (2013): they discuss how gap profiles scale differently when a planet’s Hill radius is smaller than the disk scale height. We note that our results are derived from 2D simulations, which should ideally be confirmed by 3D simulations that resolve the vertical structure of the disk to more accurately predict midplane pressure structures.

Figure 3.— Planet gap opening mass estimates based on the viscous criterion (, blue) and the Duffell (2015) criterion (, green) as a function of viscosity (top panel) and accretion rate (bottom panel). The shaded area indicates a range of a factor of 1.5 variation about the fiducial mass estimates. The red dots are FARGO simulation results based on local pressure maximum retreat. Note the Duffell (2015) criterion is a better match with the simulation results in the low and high regimes.

With our improved estimate of Vulcan planet mass set by gap opening (eq. 15), following Paper II, we eliminate the accretion rate term, , by combining it with equation 13. Previously Paper II did this using the viscous criterion gap opening mass to obtain . Now we find a revised mass versus orbital radius relation for Vulcans:


We see that this has a slightly steeper scaling of compared to the Paper II result. The implications of this and other revised IOPF predictions, such as the more general versus relation of eq. (15), for the comparison with observed STIPs will be examined in a future paper.

3.3. Formation Timescale

Now that we have improved estimates for the masses of Vulcan planets, we assess their formation timescales. Consider a protoplanetary disk with mass accretion rate of and solids to gas mass ratio of . We assume as a fiducial value. The solids are divided into two parts: (1) sub-millimeter or smaller sized dust grains that are well coupled to the gas; (2) larger, typically centimeter sized “pebbles” that can decouple from the gas flow and can be subject to substantial radial drift due to gas drag.

Pebbles are expected to arise from coagulation of dust grains, mainly via Brownian motion and turbulent mixing. The maximum steady state rate of pebble delivery to a particular radial location in the disk is if most solids are in the form of pebbles, i.e., . The accretion limited pebble supply rate is then simply:


Thus we have an estimate of the planet formation time if limited by a steady accretion rate of pebbles:


For solar type stars, estimates of the lifetimes of protoplanetary disks are to 10 Myr, with median values of 3 Myr, (e.g., Williams & Cieza, 2011; Ribas et al., 2015). Since there are usually planets in a typical STIP, this 1 Myr formation time for the first planet may give an interesting constraint on IOPF models: e.g., to form planets may take a timescale that is very similar to typical disk lifetimes. Note, the above planet formation timescale is not very sensitive to model parameters, such as opacity (), viscosity () and location ().

3.4. Pebble Supply Reservoir

The radius of the “pebble supply reservoir,” , that is needed to construct a Vulcan planet can now also be evaluated by equating this to the disk radius that encloses a solid mass that is at least as large as . We define


and show some example estimates of in Table 1, given our disk models (e.g., shown in Fig. 1). Here denotes planet formation efficiency, i.e., the fraction of solids within the reservoir that finally becomes part of planet. The size of the supply reservoir is sensitive to the value of adopted for the disk model. For the fiducial case of (i.e., potentially appropriate for dead zone conditions) and with , we find AU. However, if the value of is larger, then grows, potentially to several tens of AU, both because of the larger planet mass and the lower mass surface density of the disk.

4. Single pebble evolution model

Here we describe a model of single pebble evolution in the disk, i.e., involving radial drift due to gas drag and growth of the pebble by sweeping up small dust grains. This model is the basis of that used in the next section to predict the global evolution of the pebble population. Here we will first use the single pebble evolution model for simple estimates of the supply timescales needed to form Vulcan planets.

The pebble evolution model is based on that presented by Hu et al. (2014). It includes four different drag regimes (Epstein; and three Stokes regimes, depending on the relative radius of the pebble, to the mean free path for collisions in the gas, , and the Reynolds number, ) to evaluate the drag force . We then calculate the gas drag frictional timescale of a pebble of mass moving at speed relative to gas as :


Then the radial drift of pebble relative to gas is (Armitage, 2007):


where is power-law index of pressure radius relation in , is the dimensionless friction timescale and is the orbital angular frequency. An example of the values of as a function of pebble size and disk radius for our fiducial disk model (, ) is shown in Figure 4. We can see from this figure that efficient radial drift, i.e., when , occurs for cm-sized pebbles in the outer disk region, which will be important for the models we consider below that are based on steady injection of pebbles at this outer disk scale (see §5).

Figure 4.— Dimensionless friction timescale, , map of pebbles in our fiducial disk model (, ).
disk model (AU) (no growth) (w/ sweep-up growth)
() () (AU) () (yr, ) (yr, )
0.96 0.060 3.35, 6.87, 14.3 6.20, 7.23, 8.18 0.748, 1.17, 1.71
1.59 0.036 19.1, 40.2, 84.5 0.994, 1.22, 1.94 0.474, 0.697, 1.41
3.43 0.167 1.75, 2.95, 5.45 16.3, 24.6, 33.6 0.798, 1.35, 2.35
5.73 0.100 6.67, 13.7, 28.5 5.43, 6.37, 7.23 1.14, 1.67, 2.28
12.3 0.464 1.93, 2.90, 4.47 17.7, 29.9, 50.4 0.671, 1.14, 2.08
20.6 0.278 3.73, 6.14, 11.1 12.6, 19.8, 27.7 1.39, 2.40, 3.98

Note. – is evaluated for pebbles of initial radius of 1 mm, with the trajectory followed from to the DZIB where (for ).

Table 1Pebble Supply Reservoir Radii and Drift Timescales for Vulcan Planets

The pebble evolution model of Hu et al. (2014) involves “sweep-up growth” of pebbles by accretion of small dust grains (we add pebble-pebble coagulation in the next section). Similar models have been proposed by Safronov (1972) and also been discussed by Dullemond & Dominik (2005) in the context of dust growth during vertical settling. Windmark et al. (2012) considered such a model as a mechanism by which cm-sized “seeds” can overcome the bouncing/fragmentation barrier to form planetesimals by sweeping up smaller grains.

Like Hu et al. (2014), we assume a pebble sweeps up all the dust within its geometric cross section if it is in the Epstein drag regime, where pebble size is comparable to a gas molecule’s mean free path. However, when a pebble enters the Stokes regime, the gas behaves more like a continuous fluid, forming a pressure wake in front of the pebble. We thus assume sweep up growth stops, as small dust grains are deflected by this pressure wake and flow away following gas streamlines.

To implement this model, in each time step , the mass growth of a pebble is:


Here is mass of dust per unit volume in the disk, i.e., 1% of the gas density, , for our fiducial initial conditions. is the relative velocity between a pebble and its surrounding gas, which is given by:


We now calculate example trajectories of pebbles of initial radius of 1 mm that start at a radius in the disk equal to that of the pebble supply reservoir for a Vulcan planet, (see Figure 5). The trajectories are shown for the cases with and without sweep-up growth, i.e., the latter being for pebbles of constant radius of 1 mm. These calculations are done for six different disk models, i.e., and . In each disk, three different starting radii for are considered. These properties and the total drift times, , for all these cases are listed in Table 1.

Figure 5 and Table 1 illustrate several points. First, note the dependence of Vulcan planet masses, , and the pebble supply reservoir outer radii, , on disk properties. These planet masses range from about in the case of a low accretion rate, low viscosity disk, to about in a disk with a higher accretion rate and higher viscosity. Since disk mass surface densities are lower for higher viscosity disks, the corresponding reservoir radii also become larger in such cases. Thus varies from just 1.75 AU (actually occurring in the fiducial case if ) to almost 100 AU in the low , high , low model.

Then the drift times from these radii to the location of planet formation, i.e., , depend on whether the pebble is allowed to grow by sweep-up growth of small dust grains. Without growth, in the disks, we see that yr. It takes longer to drift inwards in the high disks, even though starting closer in, because the pebbles are in different drag regimes leading to different drag frictional timescales (eq. 21). This pattern is mirrored in the higher disks, where we see that the no growth drift times can be as short as yr even from 100 AU. Allowing pebbles to grow results in shorter drift timescales, since the pebbles approach sizes that lead to maximal drag force, i.e., minimal frictional timescales, leading to efficient radial migration, which then further enhances pebble growth. We notice that in all the six disk models considered the values of with sweep-up growth are quite similar, i.e., yr. The pebbles can grow to sizes to 10 cm.

Thus we conclude that for a wide range of disk parameters (, ), the radial drift time scale is one or two orders of magnitude shorter than , based on an estimate of a steady-state supply rate of pebbles. This reflects the extreme limit of the pebble supply rate being boosted above the accretion limited rate because of the net radial drift of pebbles with respect to gas. It also assumes that the pebbles starting from within will be able to sweep up a large fraction of the total solid mass, i.e., of dust, in this inner disk region. To advance beyond these simple estimates, we need to construct a global model of pebble evolution and supply in the disk.

Figure 5.— Radial trajectories in time of pebbles with initial radii of 0.1 cm, starting from various supply reservoir outer boundaries for Vulcan planets with (see Table 1). The black dotted lines show models with no growth of the pebbles, i.e., they keep a constant radius of 0.1 cm. The red lines show models with growth via Stokes-limited sweep-up of small dust grains. Different styles of red lines represent different ranges of pebble sizes during this growth: dashed lines cover the size range of 0.1 to 1cm; solid lines cover 1 to 10 cm. Left column shows disks with ; right column with . Top, middle and bottom rows show models with , respectively. Pebble trajectories are followed until an inner radius is reached where , , expected to be the location of the DZIB.

5. Global Pebble Evolution Model

5.1. Numerical Implementation

Here we describe the basic algorithm used in an Eulerian pebble evolution model, including radial drift and growth (sweep-up plus coagulation). We sample a discrete grid in radial distance () and pebble radius (). At any given radius, the model is intended to approximate conditions applicable to pebbles drifting inwards near the disk midplane.

We divide solids in the disk into two groups: larger particles treated as pebbles; smaller particles treated as dust. Dust represents particles smaller than a certain size threshold, for which we adopt  cm (as a fiducial choice), while pebbles are particles divided in different size bins. During radial drift, pebbles can grow by sweeping up dust (i.e., the Stokes-limited sweep-up growth model described in the previous section) and by coagulation with other pebbles.

The basic numerical approach in an “Eulerian” model is to find the fraction of material that is transported between neighboring cells within each time step. Considering pebble evolution, this fraction at radial grid cell and pebble size grid cell , is calculated as the sum of the radial drift fraction in space, and size () space . The radial drift fraction of mass moved to the next radial grid, i.e., the mass fraction of pebble size at radial grid , is the ratio between radial drift distance to radial width of the ring shaped grid:


where is the width of each radial grid, is time step, and is radial drift velocity of pebble size at the center of the grid. This “center of box” first order approximation works because the variation of disk properties within each grid is very minor, especially with a logarithmic radial grid set up.

The size evolution of pebbles is more complicated. The difference between largest and smallest pebbles within one size bin becomes larger due to the nature of sweep up growth. At size bin , the minimum size is denoted as , and the maximum is . Within one time step, a pebble with size can grow to size , by sweeping up dust of mass from eq. 23:


Similarly, a sized pebble can grow to . Between size and , there exists a size , that can grow to in this time step, . The value of is obtained from linear interpolation of and . Thus, pebbles below size will stay in current size bin, while pebbles above move to next size bin. With the assumption that pebble mass is distributed evenly within each size bin, we obtain the fractional mass in size space that grows to the next size bin:


To make the code more efficient, at each mesh in space, the time step is self adjustable, depending on the drift rate and growth rate of pebbles. The whole disk shares a larger time step for synchronization, while each mesh has its own sub time step that is a simple fraction () of the larger step. To achieve this, virtual rings are implemented between each mesh, recording fluxes from the outer mesh and supplying to inner meshes in each sub step.

Figure 6.— Example of total disk plane relative velocities between different pebbles due to radial and cylindrical drift in the fiducial disk model (, ) at 1 AU. These relative velocities are used in the pebble-pebble coagulation model.

A simple pebble-pebble interaction algorithm is also implemented. Here we ignore effects like Brownian motion, disk turbulence or vertical settling, so that only radial and cylindrical drift are considered for setting interaction velocities. When pebble size group meets group , the number of two body interactions leading to coagulation, , per unit area per time step, , is calculated as (see Birnstiel et al., 2010):


where is the coagulation cross section, is the relative velocity between pebble group j and group k, and are surface number densities of pebbles in these groups, and are pebble vertical scale heights, calculated as , and is the coagulation efficiency given by:


Thus when the relative speed is  m/s, there is a “fragmentation barrier” that prevents further pebble growth.

Figure 6 shows an example of the relative velocities occurring between different size groups at an 1 AU location in the fiducial disk with and . These range up to 40 m/s for pebbles with cm, but, for typical sizes of to 10 cm, are m/s.

The final step of modeling coagulation is that of redistributing the masses of the coagulation products. Each pebble group has a certain range of sizes (i.e., masses), and the width of this range would be enlarged by coagulation, as the product’s minimum size is set by the combination of the smallest pebbles from both groups and the maximum size by combination of the largest pebbles. So the coagulated mass will be distributed into size bins ranging from to .

The coagulation time step also follows a sub-step algorithm, in each mesh, with the condition that is chosen so that the fraction of pebbles (by number or mass) that coagulate is no larger than 0.02.

For our simulation domain, we only consider pebble evolution in the region outward of the pressure maximum of the DZIB, i.e., from about 0.1 AU to 30 AU.

The outer boundary of 30 AU is chosen as a representative outer disk scale, which contains a large enough dust reservoir to form a system of planets for the various disk models considered. However, this choice is somewhat arbitrary and, as we will see, we will mainly focus on the results of models that involve a steady injection of small pebbles at this outer boundary, so the initial reservoir of solid material is not particularly important.

The radial spatial resolutions of the model is set by 160 logarithmic divided grids, from 0.1 AU to 30 AU. We follow pebbles from a minimum particle radius of , while solids below this limit are considered to be “dust.” While  cm is the fiducial choice, we also consider models with  cm to explore the dependence of our results on this choice. Note, we use logarithmically spaced grids in pebble size space: . At the start of simulation, the dust mass surface density is set to be 1% of the gas mass surface density, . The initial pebbles are given a mass surface density of that is 1% of (in the fiducial case; variation of this parameter is also explored). The distribution of the radii of the initial/injected pebbles are set to range from to 0.3 cm with a number density versus size distribution following the equilibrium distribution of Birnstiel et al. (2015):


At the outer boundary, the disk is supplied with dust and pebbles given the steady state accretion rate and our adopted value of , and the fiducial pebble to dust mass ratio is again 1% for the injected solids.

5.2. Fiducial Model Results & Effects of Accretion Rate, Initial Pebble Distribution, & Minimum Pebble Size

Figure 7.— (a) Left: Pebble accretion rate to the inner disk (AU) region versus time for models with minimum pebble radius cm. The upper, middle and lower panels show different initial pebble to dust ratios , respectively, and the same values are also used for the injected pebble to dust ratios. Note that all models assume a constant solids (dust pebbles) to gas ratio of 0.01. In each panel, 3 different accretion rates are color coded. The solid black line shows the model of pebble evolution with only sweep-up growth for the case of the disk with . The dashed black line shows the equivalent result for the model that also includes pebble-pebble coagulation. The horizontal black dotted line shows the steady state pebble accretion rate if all solids are in the form of pebbles. Equivalent results for higher (, blue) and lower , red) accretion rate disks are shown. The lines in lighter color shading show results of disks that are initially empty of pebbles, i.e., pebbles only appear via injection at the outer boundary at 30 AU. Note the pebble accretion rates in the simulated disks can be larger than the steady state rates because of the sweep-up and delivery of the initial dust reservoir in the disk, but the asymptotic behavior at late times is towards this steady state rate. (b) Right : As (a) but now for cm.

We first present results in Figure 7 for the time evolution of the mass flux of pebbles delivered to the inner disk inner boundary () for disks with viscosity parameter . Our fiducial case adops a minimum pebble radius of cm (left column panels of Figure 7), but we also show the results for cm (right column). Our fiducial choice of initial and injected pebble to dust mass surface densities is (middle row), but we also show the results for ratios that are three times smaller (top row) and larger (bottom row). Our fiducial case has a steady accretion rate of (black/grey lines), but we explore results of varying this from (red lines) to (blue lines). We show the results of the case of pebble evolution where only sweep-up growth of dust is considered (solid lines), and then also for the case where pebble-pebble coagulation is included (dashed lines). Finally, our fiducial model involves the initial disk from 0.1 to 30 AU being populated with pebbles with the given ratio of (dark lines), however we also consider the case of an initially “empty” disk, i.e., no pebbles, only dust, with pebbles only appearing via injection at 30 AU (lighter shaded lines). The steady-state pebble supply rates assuming all solids are in pebbles are shown by the horizontal dotted lines in the panels, i.e., three lines for the three accretion rates of .

Considering the fiducial case in the left-middle panel, we see there is an initial period of pebble supply rates that are greater than the steady state rate due to delivery of solid material (mostly dust) that was initialized to be present in the disk, i.e., from 0.1 to 30 AU. Pebble-pebble coagulation leads to more efficient sweep-up of this material so that the period of elevated supply rate is shorter, but more intense. The late-time behavior of the models is asymptotic approach to the steady state supply rate (from above).

This trend is basically repeated in the higher case. However, in the lower case we see that the pebble supply rate drops below the expected steady state rate, and then slowly approaches it from below. Figure 7b shows the same models, but now for cm, which has the effect of initializing with and continuing to inject with a smaller number of larger pebbles. Similar behavior is seen, with an initial “spike” phase of elevated pebble supple rates, which are shorter, but more intense, if pebble-pebble coagulation is modeled. Now, however, the case also exhibits a pebble supply rate that falls below the steady state supply rate immediately after the spike phase.

The above behavior can be better understood, by comparing to the same results in lighter colored lines that show the models in which there are no initial pebbles in the disk: they only appear via injection at the 30 AU boundary. These models show initial elevated pebble supply rates due to delivery of the initial dust from the disk, but do not show the decremented phases where pebble accretion rate is below the steady state rate.

We conclude that the initial spike phase, due to delivery of initial disk dust, is a somewhat artificial feature of our model, dependent on the initial conditions. A subsequent decremented supply rate phase can occur in some cases if there has been very efficient sweep-up of the initial dust by an initial population of pebbles and there is a relatively long period of time needed for the injected pebbles to drift in to the inner disk. The limiting case of no initial pebbles in the disk is instructive in that it tells us the time needed for the disk to reach a quasi equilibrium state where pebble supply rate at the inner boundary equals that injected at 30 AU. For only sweep-up growth, this timescale can be Myr for the high accretion rate case, and somewhat shorter for lower cases. However, including pebble-pebble coagulation reduces these timescales. Since early high accretion rate phases are expected to last Myr, i.e., the evolution from the protostellar disk phase, this tells us that a globally self-consistent model may need to take into account the evolving disk structure due to a declining accretion rate and that a steady state may not be reached in certain circumstances.

Figure 8.— Evolution of pebble size versus orbital radius distributions from the initial condition (top row), to 0.1 Myr (2nd row), to 0.5 Myr (3rd row) to 2 Myr (bottom row). The left column shows the sweep-up only growth model. The right column shows the sweep-up plus pebble-pebble coagulation growth model. In each panel the color scale shows the distribution of the pebbles, while on the left there is a projected histogram of the sizes of those pebble inside 1 AU radii. The red solid line shows the mass surface density of pebbles. The red dot-dashed line shows the total mass surface density of solids (dust plus pebbles). When these overlap, most of the solid mass is in pebbles, which is achieved by late times in both models.

Figure 8 shows the radial profiles of the dust size distributions for our fiducial model with and without pebble-pebble coagulation. The size distribution of pebbles delivered to the inner disk is significantly larger in the latter case. The figure also shows the radial profiles of the mass surface densities of pebbles and total solids, which allows assessment of the mass fraction of solids that is in pebbles compared to dust (recall that our modeling is of the midplane region).

In models with pebble coagulation, the sweeping up of dust is not as efficient in the inner disk as in the sweep-up only model. As seen from Figure 8, the pebble mass surface density is significantly lower than the total solids mass surface density, while in the sweep-up only model, pebbles dominate the total solid mass starting even at 0.1 Myr. The growth rate of pebble coagulation scales with the square of pebble number density , while sweep up growth scales with . During the early stages in the inner disk, where the pebble density is high, fast mutual coagulation produces large, Stokes-limited pebbles without sweeping up much dust. The inner disk can maintain a high dust fraction until the pebble flux from the outer disk arrives.

Another feature worth noticing in the coagulation model is the transient behavior at 0.1 Myr. The pebbles are divided into two groups in the inner disk. There are two contributing factors: the Stokes-limited sweep-up growth forbids the further sweep-up growth of some initial pebbles in the inner disk; and the 15 m/s coagulation speed limit stops them being coagulated into larger sizes, while their mutual coagulations are very inefficient due to small cross sections and low relative velocities (see Figure 6).

5.3. Effects of Viscosity

Figure 9.— (a) Left: As Figure 7a, but now for disks with .(b) Right: As (a) but now for cm.

The appropriate value of the effective viscosity, parameterized via , is very uncertain. Here we repeat the analysis of the previous sub-section, but now for disks with . These results are shown in Figure 9. In our modeling, variation of plays a similar role as , as it mainly affects mass surface density, i.e., (see eq. 9). Higher disks thus have a smaller initial dust reservoir, i.e., in the scale of the 30 AU disk, and less efficient sweep-up growth. Thus, note that the pebble flux profile in a disk of is quite similar to the one with . Though the dust in the inner disk transfers mostly into pebbles, as in the low case, the solids in the outer disk stay mostly in the form of dust for a much longer time. This effect is visible in Figure 10, which shows the radial dependence of the pebble size distributions and the relative levels of and for the higher disks.

Figure 10.— Same as Figure 8, but showing disk with viscosity parameter . The major difference is in the outer disk: significantly more solids remain in the form of dust than in the disk with , thus resulting in a lower pebble flux after the dust in the inner disk has been swept-up.

6. Assembling a STIP

6.1. Steady Accretion Rate Disks

Figure 11.— (a) Top row: Example of STIP formation history (left and middle panels) and orbital architecture, i.e., planet mass versus orbital radius (right panel), given a pebble supply rate history for a model with cm and a steady gas accretion rate of . The results of models with only sweep-up growth of pebbles are shown in the left panel. The black solid line shows the pebble supply rate to the inner disk with , and black dashed line shows pebble supply rate with . The dot-dashed line marks the steady state pebble flux limit. The colored dots on the pebble delivery lines indicate the formation time of each planet, with blue dots for , green dots for , and purple dots for global disk but with calculation of done for , i.e., a lower value of near the DZIB. Note that in these models planet formation is assumed to start at 0.1 Myr (marked by the vertical dashed line). The middle panel shows the equivalent results, but now including pebble-pebble coagulation. (b) Middle row: As (a) but for disk accretion rates of . (c) Bottom row: As (a) but for disk accretion rates of . The disks with lower accretion rates can fail to form all five planets within the 5 Myr period considered.

We now consider the implications of these models for the timescale of assembling an entire STIP. As a specific example, we will consider the ability of disk models to form a five-planet system and the time it takes them to do so. We compare the accumulated mass from the delivered inner disk pebble flux to the required mass of the planets, i.e., set by the gap opening condition for from eq. 15.

A noticeable feature of our modeled pebble flux is the “spike” phase at the beginning of the simulations, which can transport a large amount of solids () in less than 0.5 Myr. However, this spike feature mostly reflects the adopted initial condition of pebbles and dust in the inner (AU) disk. Since these solids in reality could have already drifted in at an earlier phase of the disk’s evolution, here we ignore the pebbles delivered in the early spike phase up to 0.1 Myr. Still, we will see that our results do depend on how the initial conditions of the disks are set-up.

In the context of IOPF, building a STIP starts with the innermost “Vulcan” planet at . After that, here, for simplicity, we assume that the orbit of each subsequent planet is 20 Hill radii away from the previously formed planet, which is close to the median value from observations (Paper I) and within the range of DZIB retreat distances modeled in Paper III as being due to increased ionization via penetration of protostellar X-rays. The planet mass is calculated as using eq. 15. We simply assume all accumulated mass of pebbles will be transferred into planets, i.e., the radially drifting pebbles will be caught by the DZIB pressure trap and later accreted with 100% efficiency by the protoplanet.

The results of this modeling for our fiducial disk () are shown with solid lines in the left pair of panels in the middle row of Figure 11. Here, the first planet starts forming 0.1 Myr after the spike phase, and the fifth has completed formation within the following 0.4 Myr in the case of sweep-up only pebble evolution and within 1.5 Myr in the case including pebble-pebble coagulation (blue dots mark the time of formation of each planet). The difference is due the different pebble supply rates, with the first model more heavily influenced by the initial spike phase from the initial dust reservoir. The orbital distribution and masses of the planets are shown in the rightmost panel of the middle row (again shown as blue dots).

Equivalent models with tend to have lower pebble supply rates (dashed lines). If this value of also holds in the DZIB region, then the gap opening mass is increased, implying more massive planets. In this case only two such planets (green dots) have time to form within 5 Myr.

The top and bottom rows of Figure 11 show equivalent results for and , respectively. Higher/lower accretion rates imply higher/lower gap-opening masses, but with a sub-linear scaling. The total mass of pebbles delivered to the inner region in a given time scales more steeply with accretion rate. Thus, it becomes more difficult for lower accretion rate disk models to be able to form the entire five planet system. Overall, when comparing models with different and , we find that to form a realistic STIP, the disk should have and viscous no more than .

6.2. Evolving Disks

Lifetimes of protoplanetary disks are estimated to range from 1 to 10 Myrs (e.g., Alexander et al., 2014; Espaillat et al., 2014). The effects of an evolving disk have been previously considered in dust evolution models. For example, Birnstiel et al. (2010) set up a disk model subject to viscous evolution and gravitational instability, starting from the very early infall phase. It gives a relatively constant accretion rate starting from 0.15 Myr to 1 Myr. Other methods involve adopting simple functional forms for the evolution of the accretion rate, such a power law or exponential decay (e.g., Hueso & Guillot, 2005; Bitsch et al., 2015). The disk modeled by Bitsch et al. (2015) has an accretion rate that evolves as:


with and yr. We adopt the same functional form, but with and yr, so that the accretion rate declines by a factor of 10 in 1 Myr.

In our numerical modeling, we update the disk properties, i.e., the accretion rate, every 1000 years. Note, in such models of exponentially declining accretion rate, the disk evolution is not fully self consistent for its accretion history, i.e., the declining does not match the viscous radial drift of gas. Thus these models should be viewed as being simple “toy” models of disk evolution. Note also, when we force the gas disk properties to evolve with an exponentially declining rate, we have a choice about how we treat the associated dust that may be mixed together with the gas. We consider two extreme cases: Case A: the dust content of the disk, i.e., , does not evolve with the declining gas, i.e., it would remain at a steady mass surface density, except for the effects of radial drift (due to the assumed viscosity) and sweep-up by pebbles; Case B, the dust content of the disk, i.e., , is decreased with the same exponential decay rate of the gas, with this decay factor applied multiplicatively on top of any other changes due to radial drift and sweep-up by pebbles.

In principle, there are now several parameters that may be varied in these models: , , (i.e., the time when the first planet starts to form), in addition to other model parameters already introduced, such as and . Furthermore, there are the choices of: Case A or B, introduced above; whether or not pebble-pebble coagulation is included in the modeling; and whether or not the disk model is initially “empty” of pebbles. It is not our intention here to carry out a full exploration of the parameter space of all these models. Rather, here we present several example models of representative cases, in particularly focussing on the examples of parameter combinations that produce realistic STIPs within reasonable timescales.

Figure 12.— (a) Top row: Example of STIP formation history (left and middle panel) and orbital architecture (right panel) in an evolving disk (treated as Case A, see text), using same color coding and styles as Figure 11. Here, planet formation also starts at Myr. The only addition is the lighter shaded points in the right panel represent planets formed in the “with coagulation” model. (b) Bottom row: As (a), but now for Case B, i.e., dust content of disk is also scaled with the exponential decay function.
Figure 13.— (a) Top row: Same as Figure 12a, i.e., Case A models, but now with planet formation starting at Myr. (b) Bottom row: As (a), but now for Case B, i.e., dust content of disk is also scaled with the exponential decay function.

Following the format of Figure 11, Figure 12 now shows example models with , yr,