Save the Planet, Feed the Star: How Super-Earths Survive Migration and Drive Disk Accretion
Jeffrey Fung Eugene Chiang



Department of Astronomy, University of California at Berkeley, Campbell Hall, Berkeley, CA 94720-3411
1NASA Sagan Fellow

email: jeffrey.fung@berkeley.edu



Two longstanding problems in planet formation
include (1) understanding how planets survive migration,
and (2) articulating the process
by which protoplanetary disks disperse—and in particular
how they accrete onto their central stars.
We can go a long way toward solving both problems 
if the disk gas surrounding planets
has no intrinsic diffusivity (“viscosity”).
In inviscid, laminar disks, a planet readily repels gas
away from its orbit.
On short timescales,
zero viscosity gas accumulates
inside a planet's orbit to slow Type I migration by orders 
of magnitude. On longer timescales,
multiple super-Earths (distributed between,
say, ∼0.1–10 AU) can torque 
inviscid gas out of interplanetary space, either inward 
to feed their stars, or outward to be blown away in a wind.
We explore this picture with 2D hydrodynamics simulations
of Earths and super-Earths embedded in inviscid disks,
confirming their slow/stalled migration even under
gas-rich conditions,
and showing that disk transport rates
range up
to ∼10^-7 and scale as 
Ṁ∝Σ M_ p^3/2, 
where Σ is the disk surface density 
and M_ p is the planet mass.
Gas initially sandwiched between two
planets is torqued past both into the inner
and outer disks.
In sum,
sufficiently compact systems of super-Earths can
clear their natal disk gas, 
in a dispersal history
that may be complicated and non-steady, but which conceivably
leads over Myr timescales to large gas depletions 
similar to those characterizing transition disks.








§ INTRODUCTION




Protoplanetary disks have two jobs: make planets and 
feed their host stars. The first task is frustrated
by migration: disk torques force 
planetary orbits to decay <cit.>,
evacuating the very regions where planets
are observed in abundance 
<cit.>.
Fulfilling the second task requires a mechanism
to transport away the disk's angular momentum.
Magnetic torques are promising but depend
on seed fields of uncertain provenance <cit.>.

<cit.> proposed that the two problems are actually
one: that planets themselves—if they can survive migration—can 
provide an effective source
of disk viscosity by exciting density waves that transport
angular momentum outward. <cit.> emphasized that
such planets must be massive enough to open gaps and
avoid Type I migration.
They focused on giant Jupiter-mass planets, a demographic
that is now understood to be rare <cit.>.
In this paper, we turn our attention to super-Earths: bodies
of mass 1–10 M_⊕ that have been
discovered by Kepler to be
relatively commonplace <cit.>.

Can super-Earths avoid Type I migration? Yes—if
their disks are sufficiently inviscid. The
dependence of the Type I drift rate on disk viscosity
is perhaps under-appreciated, as it is not explicitly
called out in the typically quoted Type I formula
(see, e.g., ).
Crucially, without an intrinsic disk viscosity to 
smooth away the planet's perturbations
to the disk's surface density,
a pile-up of disk material ahead of a migrating planet
exerts a “feedback” torque that
slows and can even stall migration
<cit.>.
<cit.> calculated
that for inviscid
disks in which planet-driven waves dissipate 
by steepening into shocks, the critical
planet mass above which Type I migration shuts off is:

    M_ cr≃
    4 (h_ p/r_ p/0.035)^3(M_∗/M_⊙)
    (Σ_ p r_ p^2 / M_∗/10^-3)^5/13 M_⊕  ,

where r and h are the disk radius and scale height,
Σ is the disk gas surface density, M_∗ is the
central stellar mass, and the subscript p
indicates evaluation near the planet's position. 
<cit.> and <cit.> have performed numerical
simulations supporting the analytic calculations by 
<cit.>, and confirming that super-Earths
in low-viscosity disks migrate much more slowly (and erratically)
than is predicted by Type I.

A planet of mass M_ p
at r_ p drives a disk
mass transport
rate Ṁ at
distance r of

    Ṁ(r) = -2F_0/l r ∂φ(r)/∂ r

where

    F_0 = Σ_ p r_ p^2 l_ pΩ_ p(M_ p/M_*)^2 (h_ p/r_ p)^-3

measures the total angular momentum carried away per time
by planet-driven waves (a.k.a. the total integrated
one-sided Lindblad torque),
Ω is the orbital frequency, and
l = Ω r^2 is the specific angular momentum.
The dimensionless function φ(r) 
describes how waves, as they travel away from the planet,
damp with distance, depositing their
angular momentum to disk gas 
and thereby propelling material
radially.
From <cit.>,

    φ∼( M_ p/M_ thermal)^-1/2( |r_ p - r|/h_ p)^-5/4 ,

valid for M_ p≲ M_ thermal≡ (h_ p/r_ p)^3 M_∗ and
φ≲ 1 (i.e., distances far enough from the planet
that the waves are dissipating in weak shocks), and where
we have ignored order-unity constants and
all radial variations in h, Σ, and gas sound speed.
It follows that

    Ṁ(r) ∼ sign(r-r_ p) Σ_ p r_ p^2 Ω_ p(M_ p/M_∗)^3/2(h_ p/r_ p)^-5/2(|r_ p-r|/h_ p)^-9/4  .

At r < r_ p, Ṁ as given by
Equation (<ref>) is negative
(mass flows inward), and vice versa; a planet
tends to repel material away from itself.
To avoid minus signs, we will ignore this formal
sign convention so that all our reported values for Ṁ
will be positive and understood to be inward unless
otherwise indicated.

Note how Ṁ∝ M_ p^3/2 and not M_ p^2.
Although the total Lindblad
torque scales as M_ p^2 (Equation <ref>),
that torque is distributed over a distance that increases
with decreasing M_ p (as M_ p^-2/5, as can be
seen by solving for |r_ p - r| in terms of M_ p
at fixed φ in Equation <ref>). 
Thus at fixed distance away from the planet, Ṁ increases
with M_ p with a power less than 2.

Inserting M_ p  = 10M_⊕
and other nominal parameters (for r < r_ p) into
(<ref>) yields

    Ṁ(r) ∼ 10^-8   (Σ_ p r_ p^2/10^-3 M_⊙) (2π/Ω_ p/1  yr)^-1(M_ p/M_*/3 × 10^-5)^3/2(h_ p/r_ p/0.035)^-1/4
       ((r_ p - r)/r_ p/0.5)^-9/4  ,

comparable to accretion rates measured for classical T Tauri stars
<cit.>.
Note how weakly Ṁ depends on h_ p/r_ p,
underscoring how Ṁ does not scale simply
as the total Lindblad torque in eqn:torque 
(which scales as (h_ p/r_ p)^-3),
but depends also on the distance over which that torque
is exerted, as we have described above.

The above considerations indicate that with super-Earths
we might have our cake (survive migration) and eat it, 
too (drive disk accretion). Of course, a single 
super-Earth is insufficient because its reach is too 
short (Ṁ drops as |r-r_ p|^-9/4).
Multiple super-Earths are needed to shuttle the
accretion flow from distances of a few AU
down to the stellar radius. Reality will be non-steady
and likely messy (see, e.g., Figure 4 of
), with material between
adjacent planets having a fate that is not obvious:
does the sandwiched gas drain inward,
or does the inner planet hold back material
pushed inward by the outer planet?
And to what extent do super-Earths migrate with
the accretion flow they drive?

Here we explore these questions
using fully non-linear, 2D
hydrodynamical simulations of super-Earths embedded
in inviscid disks. We measure the migration
histories r_ p(t) and accretion rates Ṁ
in simulations containing 1 or 2 super-Earths,
experimenting with varying the disk surface density and
the planet mass to test eqn:ana_mdot_wnum.
Our numerical methods are given in Section <ref>.
Results are presented in Section <ref> and
placed into broader context in Section <ref>.




§ NUMERICAL METHOD




We use the graphics processing unit (GPU) accelerated
hydrodynamics code  <cit.> to perform
2D simulations of disk-planet interactions.
It is a Lagrangian-remap shock-capturing code that uses the
piecewise parabolic method <cit.> to solve the
continuity and momentum equations:

    DΣ/ Dt   = -Σ(∇·𝐯)  ,
     D𝐯/ Dt   = -1/Σ∇ p + 1/Σ∇·𝕋 - ∇Φ ,

where Σ is the gas surface density, 𝐯 the
velocity field, p the vertically averaged gas pressure,
𝕋 the Newtonian stress tensor, and Φ the
combined gravitational potential of the star and the planet(s).
We use a globally isothermal equation of state:
p=^2Σ with a spatially constant sound speed 
=0.035 v_ K, 1 AU≃ 1  km s^-1
where v_ K, 1 AU is the Keplerian velocity
at 1 AU around a 1 M_⊙ star. (This  corresponds
to a disk temperature of 300 K assuming a mean molecular weight
of 2.34.)

In a polar coordinate system (radius r, azimuth ϕ)
centered on the star,

    Φ   = -GM_*/r + ∑_i=1^N_ pΦ_ p,i
    Φ_ p,i   = -GM_ p, i/√(r^2 + r_ p, i^2 - 2rr_ p, icosϕ_i' + r_ s, i^2) + GM_ p, i rcosϕ_i'/r_ p, i^2

where G is the gravitational constant,
M_∗ = 1 M_⊙ is the stellar mass,
the subscript i labels each planet,
N_ p is the
total number of planets, 
M_ p is the planet mass, 
Φ_ p the planet's gravitational potential,
r_ p the planet's radial coordinate, 
 the smoothing length
of the planet's potential, and ϕ' = ϕ-ϕ_ p the
azimuthal separation from the planet.
The stress tensor 𝕋 is proportional to the
kinematic viscosity ν. Most of our simulations are of 
inviscid disks with ν=0. For our viscous disk
simulations, we use ν=α h, where the Shakura-Sunyaev
parameter α = 0.001, h = /Ω_ K is the
local scale height, and Ω_ K = √(GM_∗/r^3)
is the Keplerian orbital angular
frequency. At r = 1 AU, h/r = 0.035.
We set = 0.5 h,
as is appropriate for 2D simulations <cit.>.

A given planet feels the gravitational force from the star,
the disk, and other planets.
The disk force on the planet is calculated by direct summation
over all mass
elements in the disk, with the “background”
axisymmetric component of
the disk surface density subtracted off. Because the disk does not
feel its own gravity at all
(i.e., we ignore disk self-gravity;
see equation <ref>), eliminating
this axisymmetric component in the disk-planet forcing
improves consistency between the motions
of the planets and the disk. Planet migration should be minimally
affected by this procedure, since the background
component of Σ exerts
no torque. 
Spurious forces arising from within the planet's Hill
sphere
are
sometimes a concern if
this
region is under-resolved.
The Hill radius, = (M_ p/3M_*)^1/3, 
ranges from 0.3 to
0.6 h_ p, similar to the smoothing length .
We have
verified that the
torque generated 
within a radius of
0.5  from the planet is negligible, and so
we do not excise
the Hill sphere
in force calculations.
The planets' motions are integrated using a
kick-drift-kick leapfrog scheme, with the drift step occurring
synchronously with the hydrodynamics step; i.e., the planets'
positions are linear in time within a hydrodynamics step.






 §.§ Initial and boundary conditions, and grid parameters



ccccccccc
 Model Parameters
Model #   M_ p (M_⊕)   r_ p,1 (AU)   r_ p,2 (AU)   Σ_0 (g cm^-2)     α     r_ in (AU)    t_ end (years)
1     10     1        –       8.5× 10^3     10^-3     0.4     700

2     10     1        –       8.5× 10^3     0             0.4     5000

3     10     1        1.2      8.5× 10^3     0             0.4     5000

4     10     0.75     1.05     8.5                    0             0.3     2000

5     3      0.75     1.05     8.5                    0             0.3     2000

6     1      0.75     1.05     8.5                    0             0.3     2000
t_ end is the end time of a simulation, in units where the Keplerian orbital period at 1 AU is 1 year. Also, r_ p,1 and r_ p,2 are merely the initial planet locations at t=0; the planets are completely free to migrate in the simulations.


tab:setup lists the parameters used by our
6 models.
The disk is initialized with a power-law surface density:

    Σ =  Σ_0 (r/ AU)^-3/2  .

We consider both gas-rich disks having
Σ_0=8.5×10^3  g cm^-2
resembling the minimum-mass
extrasolar nebula (),
and gas-poor disks having a surface density
1000× lower. The initial velocity field
is axisymmetric and Keplerian,
with corrections from gas pressure:

    Ω = √(Ω_ K^2 + 1/rΣ d p/ d r) .

One planet, whose mass is increased gradually over the
first 10 yr of the simulation
to the full value of M_ p (either
1, 3, or 10 M_⊕), 
is placed initially at
r = r_ p, 1 (either
1 or 0.75 AU) 
and ϕ_ p, 1 = π.
In two-planet models, we place a second planet
of equal mass to the first at
r = r_ p, 2 (either 1.2 or 1.05 AU)
and ϕ_ p, 2 = π initially.

Our simulation grid spans the full 2π
in azimuth, and extends from an outer radius
of 1.8 AU to an inner radius r_ in
that equals either 0.4 or 0.3 AU depending
on whether r_ p,1 = 1 AU or 0.75 AU
(see tab:setup). 
Grid dimensions are 800 (r) × 3200 (ϕ) when
r_ in = 0.4 AU,
and 960 × 3200 when r_ in = 0.3 AU.
Cells are spaced logarithmically in radius and
uniformly in azimuth.
Our choices yield a resolution of ∼18 cells
per scale height h
in both directions at r = 1 AU
(similar to the resolution of ).
Simulations at twice our standard resolution
did not produce
significant changes in either planet migration
or disk accretion rate for the first 100 yr.
We also tested our inviscid disk model without a planet,
and found that the numerical noise in |Ṁ| 
was about 3 orders of magnitude below planet-driven
disk accretion rates,
corresponding to a numerical viscosity of α<10^-5.

Radial boundary conditions require special care in this study.
After experimenting with a few ways to measure disk accretion rates,
we found that the most stable method was to track the total
disk mass within a cylinder of radius 0.6 AU—a distance
intermediate between
the innermost planet and the inner disk boundary—while
preventing mass from leaving the grid.
We adopt “zero flux” boundary conditions where mass and momentum
fluxes across the inner and outer disk edges are always zero.
In , this is achieved by solving a special 
Riemann problem at the boundaries, one where no wave
travels toward the simulation domain, and where the radial velocity
outside the domain is always zero. This implementation
conserved the total mass within the simulation domain
to numerical accuracy.
The accretion rate Ṁ at r = 0.6 AU
is calculated by following over time
the disk mass enclosed, M_ 0.6  AU(t).
Because the function M_ 0.6 AU(t) fluctuates strongly, we fit
independent lines to segments of data each lasting 20 yr,
taking Ṁ from the best-fitting slopes.

As a planet repels material away from its orbit, our boundary
conditions result in gas piling up at the inner and outer boundaries.
Our results can only be trusted to the extent that
these boundary pile-ups do not interfere with planet migration and disk
accretion. We therefore limit ourselves to studying only the
first few thousand years of planet-disk interactions, before 
boundary effects become too large.




§ RESULTS




We assess to what extent planets migrate in inviscid disks
(<ref>),
and study how planet-driven accretion rates evolve with time and
depend on disk and planet masses (<ref>).
For planet migration, models #1–3 demonstrate differences
between viscous and inviscid disks, and between single-planet and two-planet
systems. For planet-driven accretion, we vary disk and planet masses
in models #3–6 to test eqn:ana_mdot_wnum.



 §.§ Planet migration











fig:mig plots the orbital evolution of 
planets, each of mass 10 M_⊕, 
in the gas-rich disk models (#1–3). Overplotted
for comparison is the trajectory expected from integrating
the Type I migration rate,

    ṙ_ p, Type  I = -2 C r_ pΩ_ p(Σ_ p r_ p^2/M_ p) (M_ p/M_*)^2 (h_ p/r_ p)^-2  ,

using the unperturbed surface density law in
eqn:initial_sigma to evaluate Σ_ p
(using the actual surface density in the viscous
disk simulation #1 would give practically identical
results, since gaps do not form in that model).
Three-dimensional simulations suggest that
C∼ 2–3 for our given disk
profile <cit.>; the blue
dashed curve in fig:mig uses C=2.
fig:sigma displays
azimuthally averaged surface density profiles at various
epochs in gas-rich, inviscid disk models #2 and #3,
with planet locations marked.

In agreement with the simulations of <cit.>
and <cit.>, the planet
migration rate in our viscous disk (model #1)
is similar to the Type I rate (punctuated by what appear
to be episodes of even faster Type III migration;
), and
much slower in inviscid disks (models #2
and #3). The initially rapid migration 
seen in the inviscid simulations at t ≲ 200 yr
is a transient that decays after disk surface densities
adjust to planetary Lindblad torques, i.e., 
after the surface density pile-up ahead of the planet 
attains a fractional amplitude on the order of unity.
After this initial adjustment period,
migration slows and even stalls at times, with radial positions
changing by ∼10–30%, or less, over kyr timescales.

In the single-planet, gas-rich simulation (black curve
in fig:mig,
model #2), the planet journeys slowly inward for the first
∼2000 yr and practically stops in the mean
from t ≃ 2000–3500 yr,
as disk gas that the planet has pushed inward to
r ≃ 0.7–0.8 AU piles up (blue and
green curves in the left panel of fig:sigma)
and stymies further migration.
At t ≃ 3500 yr, the planet experiences a sudden
drop in orbital radius; we traced this drop to
a close encounter between the planet and a vortex formed at its outer gap
edge at r ≃ 1.1 AU (see fig:vortex).
Some time after the encounter, the vortex gradually disperses,
completely decaying away by the end of our
simulation at 5000 yr (fig:endpoint, left panel).
Similar planet-vortex interactions were found
in simulations by <cit.> and <cit.>.
Thereafter, the planet's migration returns to its
near-zero mean pace.

In our two-planet, gas-rich simulation (#3),
the outer planet stalls for the first ∼3000 yr (magenta
curve in fig:mig), apparently trapped
at a local surface density maximum created
by the inner planet whose forcing 
dominates: see the magenta curve
in the right panel of fig:sigma,
and note how similar it is to the corresponding
magenta curve in the left panel for the single-planet case.
Gas pushed inward by the outer planet strengthens the
torque on the inner planet
and forces the latter
to migrate inward by ∼20% over the same time period
(red curve in fig:mig);
contrast this behavior with the stalling observed
in the single-planet case (black curve). 
Eventually, at t ∼ 4000 yr,
the outer planet disperses the surface density
maximum in its vicinity, and proceeds to migrate
slowly inward, slowing down near t∼ 5000 yr
as it runs into sandwiched gas. Meanwhile,
the inner planet ultimately comes to a near halt
in much the same way that it does in the single-planet 
simulation, having run into material that has piled up
just interior to its orbit.
The right panel of fig:endpoint
shows the final surface density distribution.
We emphasize that this pile-up is physical as it
is located at r ≃ 0.6 AU, away
from the inner grid boundary of the simulation at
r_ in = 0.4 AU.
The latter location has
its own separate pile-up, which does not grow
to significance over the limited duration of our simulations.
The same statement applies to the outer grid boundary.

We note that in none of the two-planet simulations did we observe
the formation of a vortex like the one seen in our single-planet
simulation. This difference might be physical,
and deserves attention in future studies of planet-vortex 
interactions.

All other factors being equal, lower disk masses should lead
to even slower planetary migration rates. 
This is confirmed in our gas-poor simulations (models #4–6)
which exhibit no measurable change in planet mean
radial positions. Thus the gas-poor simulations
can be used to diagnose disk accretion rates
without the complicating effects of planetary migration,
as we discuss in the next subsection.




 §.§ Disk accretion







fig:acc shows disk accretion rates
as functions of time for all our inviscid models. 
As described in <ref>, the accretion
rate Ṁ is measured by tracking the build up
of disk mass inside r = 0.6 AU, a location
interior to the planets at all times. 
Later, in <ref>,
we track the movement of mass initially
between two planets.

We begin by checking whether
our simulations are compatible with the analytic expectation
for Ṁ given by eqn:ana_mdot_wnum.
The comparison is best made at early times of the simulation,
t ≲ 1000 yr,
before radial surface density profiles
become too distorted.
For the single-planet, gas-rich
model #2, eqn:ana_mdot_wnum yields
Ṁ∼ 3 × 10^-8. This
prediction is within a factor of ∼2 of the 
simulated result at early times
(left panel of fig:acc).
As for the corresponding two-planet model #3,
we expect from (<ref>)
that Ṁ should be only fractionally
larger than for the single-planet case; the
second planet is farther removed from where
we measure Ṁ (r = 0.6 AU),
and so makes only a ∼50%
contribution to the accretion flow there as
compared to the inner planet.
This is approximately consistent with fig:acc.

For our gas-poor model #4, eqn:ana_mdot_wnum predicts
Ṁ∼ 1.5 × 10^-10, again
within a factor of 2 of the simulated result (right panel of
fig:acc). Scaling the planet mass M_ p
down by a factor
of 10 from models #4 to #6 should, according to
eqn:ana_mdot_wnum, reduce Ṁ by a factor of
10^3/2≃ 30. By comparison, fig:acc
shows a factor of ∼20
decrease between these two models; we consider
this acceptable agreement with the analytic expectation.
In summary, our simulations
support the various functional dependencies
predicted by eqn:ana_mdot_wnum
to within a factor of 2.

At later times, t ≳ 1000 yr, 
we observe time variability in Ṁ caused by
the deepening of planetary gaps, and by
planet migration. These variations are limited
to factors of a few. 
The simulations easiest to interpret are models #4–6 (right
panel of Figure <ref>) which have too little disk gas
to drive planet migration. The initial gradual decline in Ṁ
seen in model #4 is caused by the deepening
of gaps opened by its 10-M_⊕ planets;
over the course of 2000 yr, the gas density in the immediate vicinity
of the planets decreases by a factor of ∼5 for the inner planet 
and by a factor of ∼3 for the outer one.
Models #5 and #6 exhibit steadier accretion rates,
as their planets have lower masses which are less effective
at opening gaps.

More complicated behavior is seen in the gas-rich simulations
where planets migrate more appreciably.
Comparison of Figures <ref> and <ref>
reveals that increases in Ṁ can be traced to planets
moving inward, either gradually, as in the first 2000 yr
of models #2 and #3, or suddenly, as in the
planet-vortex encounter at t ≃ 3300 yr in model #2.
Decreases in Ṁ correspond to planets opening gaps
upon moving to new locations.



  §.§.§ The fate of gas initially sandwiched between planets


Disk accretion driven by planets would be impractical
if material residing between planets were unable to escape.
We track this sandwiched gas in simulations
#4–6, each containing a pair of planets which migrate
negligibly. We assign each gas parcel a “passive scalar” η
that equals 1 for gas initially located between 
r = 0.8 and 1 AU (between the two planets), and is 0 everywhere else.
Gas elements carry η as a conservative quantity.

fig:eta follows the η-tagged gas by plotting 

    ∂M_η/∂ r = ∫^2π_0ηΣ r dϕ  ,

vs. r at various times.
We find that the sandwiched gas 
is torqued both inward and outward,
escaping in roughly equal amounts
to the inside of the inner planet
and to the outside of the outer planet.
The opposing torques from the
two planets do not in general cancel.
The opposing torques from the  
two planets do not in general cancel, although 
in model #4 containing the highest mass planets,
some gas does concentrate along the midline between
the planets in a “shepherded” ring, resulting in less
mass leaking out of the sandwiched region. 
Looking at models #4–6 in Figure <ref>,
we see no clear trend between the rate at which sandwiched gas escapes
and planet mass. There seems to be a complicated confluence of effects
in the sandwiched region. Lindblad torques act to shepherd some gas
while also opening gaps that reduce the local gas density; and
co-orbital torques allow gas to escape via horseshoe orbits,
whose libration times scale only weakly with planet mass  
(t_ lib∝ M_ p^-1/2; ).
A detailed analysis is deferred to another paper; for now,
we conclude that, at least for comparable mass planets
with orbital spacings like the one we have assumed,
the inner planet presents a porous barrier
to material pushed inward by the outer planet.
Apparently gas that is pushed by the outer planet
toward the inner planet can be shuttled
past the latter on horseshoe orbits (and vice versa).







§ SUMMARY AND DISCUSSION




Using hydrodynamical simulations, we have demonstrated
that super-Earths in inviscid disks can simultaneously avoid type I migration
(fig:mig) and promote disk accretion (fig:acc)
by driving density waves. Disk accretion rates
measured from our simulations verify analytic predictions
(eqn:ana_mdot_wnum) to within a factor of 2 .
We observed gap opening and planet migration in inviscid disks
to be modest and to introduce order-unity effects on the disk
accretion rate. We also found in our two-planet simulations
that material initially sandwiched between two planets leaks past both
into the innermost and outermost disks (fig:eta).

Our models omit a number of effects. Many of these are not overly
concerning. Although our simulations are 2D, no substantive difference
between 2D and 3D treatments of planet-disk interactions 
in viscous disks has been reported
vis-à-vis gap opening <cit.>
or planetary torques <cit.>. Our neglect of disk self-gravity
should be an excellent approximation, as
the Toomre Q-values of our disks greatly exceed unity.
Our planets are not allowed to accrete gas,
but super-Earths/sub-Neptunes
are inferred observationally to have only modest amounts
of gas—less than 10% by mass—acquired gradually
over the entire disk lifetime <cit.>.

More interesting frontiers to pursue include
incorporating disk thermodynamics, 
as radiative cooling and differential heating across gap walls
are thought to materially affect planet-disk interactions
(e.g., ; ).
Of course, extending the durations of the simulations,
and including more planets with different orbital
architectures, would also be welcome, for greater realism and
to enable more direct connections with observations.
Closer study of gap depths is warranted;
we observed surface density contrasts only on the order
of unity (fig:sigma), in clear deviation from scaling relations
derived from viscous disks <cit.>,
and surprising insofar as less viscous gas should be
less effective at diffusively back-filling gaps.
The shallowness of the gaps is due partly to the planets migrating 
and re-starting the gap-opening process at each new radial location <cit.>.
But how much of it is due to the limited duration of our
simulations (≤ 5000 yr), or to hydrodynamical instabilities
like the Rayleigh instability <cit.>, remains
to be worked out.
Finally, survival against planetary migration
is not guaranteed: planets with mass M_ p > M_ cr
(eqn:m_cr) stall but less massive planets do not.
The question is whether rocky planets can
coagulate fast enough to cross the M_ cr threshold
before they succumb to migration.

Our results support the proposal
by <cit.> that predominantly rocky
planets—super-Earths and Earths—can solve, or at least help
to solve, the problem of how protoplanetary gas disks ultimately disperse.
Given that a single planet can push gas over a lengthscale of
approximately half its orbital radius,
shuttling gas from 5 AU down to 0.1 AU would require
about 6 super-Earths distributed in roughly equal logarithmic intervals
across this distance. Such planet multiplicities
are reasonable, given the profusion of super-Earths/sub-Neptunes
discovered by Kepler <cit.>.
To be sure, a disk accretion flow driven by a planetary system
will be unsteady, changing not only on secular, Myr-long timescales,
but also on much shorter ones, with mass alternately
accumulating and dispersing
in interplanetary space, as we have seen in our simulations
(fig:sigma).
No matter how complicated the history, however,
all gas must ultimately be torqued out of
sufficiently compact planetary systems.
It may be torqued by planets
so far inward that turbulence driven by the
magneto-rotational instability, activated in the innermost
regions which are sufficiently thermally ionized
<cit.>,
takes over the job of disk accretion onto the host star.
Or it may be torqued by planets
so far outward that it escapes from the system
altogether in a photoionized wind <cit.>.

How massive a gas disk can a set of super-Earths drain?
A first consideration is that the disk surface density
can not be so large that the embedded planet
mass M_ p < M_ cr∝Σ_ p^5/13,
lest the planet migrate away. Based on the typical
parameters listed in eqn:m_cr, a disk
containing Σ_ pr_ p^2 ∼ 10^-3 M_∗∼ 300 M_⊕
of gas could be evacuated by ∼6 super-Earths
weighing a total of ∼30 M_⊕.
In such an initially gas-rich environment, we anticipate
the planets would migrate to and fro by a few tens of percent
in orbital distance (fig:mig).
As the disk drains in the long term, whatever
slow and erratic migration the planets
undergo diminishes.
Most of the angular momentum
of interplanetary gas would be transported
to the outermost disk, exterior to all the planets, 
either by Lindblad torques or by direct advection.

The large cavities of transitional disks may have been
excavated over time by families of super-Earths. 
Observationally, gas densities inside cavities
can be suppressed relative to their values
outside by two to four orders of magnitude
(e.g., ; ).
To reproduce these strong depletions,
appeal is commonly made to giant planets in viscous 
α-disks that can open deep gaps <cit.>.
But an alternative interpretation is that the cavities
have been eroded gradually over time by much smaller mass planets
in inviscid disks. So far as we have measured in our simulations,
such planets open gaps having only order-unity
surface density contrasts. Nevertheless, given
sufficient time, they can drain interplanetary gas
by orders of magnitude.
Because planet-driven accretion rates Ṁ scale linearly
with gas surface densities Σ, we have Σ̇∝ -Σ
which implies exponential decay of the gas content.
If it takes t_ e-fold∼ 10^5 yr
to reduce a total disk mass of 10^-3 M_⊙ by a factor of e
(this assumes a contemporaneous mass transport rate
of 10^-8),[This transport can be inward
or outward—it should not matter as long as the region occupied
by the planets is monotonically drained of gas over time.]
then it takes 7 t_ e-fold∼ 7 × 10^5 yr to reduce
it by a factor of 1000.
Another feature of this picture is that if 
super-Earths drain disk mass faster than it takes for their
nascent atmospheres to
cool and acquire more mass,
they may be able to 
forestall
runaway accretion <cit.>.
Depleting the local disk density by a factor of ≳ 100
(relative to the minimum-mass extrasolar nebula) 
over timescales of ∼1 Myr suffices to keep
the gas-to-solids mass fraction of super-Earths
≲ 10%, in accord with observations <cit.>.

One potential problem with this scenario
is that it predicts mass accretion rates to lower in proportion
to disk gas densities. Although some transitional
disks do have low accretion rates <cit.>,
others do not, with a few having Ṁ as high
as 10^-7 (e.g., ; ;
). Even so, system-to-system
variations in orbital architectures,
particularly in the masses of the orbiting
companions, could help to resolve
this problem. We have focused here on super-Earths because
they are commonplace, but in principle gas giants and perhaps brown
dwarfs or even low-mass stars may serve in their stead.[The 
case of the transitional disk HD 142527 is especially intriguing:
its host star accretes at a rate of Ṁ∼ 10^-7
<cit.>,
and its cavity contains a 0.2 M_⊙ companion
highly inclined to the disk <cit.>.]

That disk gas is intrinsically inviscid (laminar)
is suggested on other, independent grounds.
Large-scale asymmetries in transitional disks (e.g.,
; ; )
have been interpreted as vortices (e.g., ;
), but these vortices
are spawned only in low-viscosity disks. 
<cit.> found that the Shakura-Sunyaev viscosity
parameter α needed to be 10^-4 or lower before vortices could
grow from sharp density gradients.
In other news, attempts to detect turbulence in the outer portions
of disks using molecular line observations 
have so far come up empty-handed (;
Flaherty et al. 2017, in preparation).
And <cit.>, in a systematic analysis of disk
accretion rates and masses, suggests that accretion
may not proceed viscously (i.e., diffusively),
but may be enabled instead by 
spiral density waves and/or disk winds.
All these recent developments, in addition to our present work, suggest 
that the reason the community has not discovered
a robust explanation for a non-zero α in protoplanetary
disks is that none exists: that in fact such disks are for the most part
inviscid, and accrete 
primarily by the action of gravitational torques,
exerted either by disk gas itself at early times
(e.g., ; ),[In 
self-gravitating disks, characterizing transport in terms
of a non-zero α is commonly done, but mostly
for convenience. Gravity is a long-range force and
not naturally captured within a local theory
like the one defining α.]
or by planets at late times.

We thank Pawel Artymowicz, Ruobing Dong, Eve Lee, Hui Li, Zhi-Yun Li, Frédéric Masset, Norm Murray, Ruth Murray-Clay, Sijme-Jan Paardekooper, Roman Rafikov, Yanqin Wu, Zhaohuan Zhu, and an anonymous referee for encouraging discussions and helpful feedback.
This work was performed under contract with the Jet Propulsion
Laboratory (JPL) funded by NASA through the Sagan Fellowship Program
executed by the NASA Exoplanet Science Institute.
EC is grateful for financial support from NASA and NSF.


apj



