Issue |
A&A
Volume 497, Number 3, April III 2009
|
|
---|---|---|
Page(s) | 869 - 888 | |
Section | Planets and planetary systems | |
DOI | https://doi.org/10.1051/0004-6361/200811265 | |
Published online | 05 March 2009 |
Planet formation bursts at the borders of the dead zone in 2D numerical simulations of circumstellar disks
W. Lyra1 - A. Johansen2 - A. Zsom3 - H. Klahr3 - N. Piskunov1
1 - Department of Physics and Astronomy, Uppsala Astronomical Observatory, Box 515, 751 20 Uppsala, Sweden
2 - Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
3 - Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
Received 31 October 2008 / Accepted 23 February 2009
Abstract
Context. As accretion in protoplanetary disks is enabled by turbulent viscosity, the border between active and inactive (dead) zones constitutes a location where there is an abrupt change in the accretion flow. The gas accumulation that ensues triggers the Rossby wave instability, which in turn saturates into anticyclonic vortices. It has been suggested that the trapping of solids within them leads to a burst of planet formation on very short timescales.
Aims. We study in the formation and evolution of the vortices in greater detail, focusing on the implications for the dynamics of embedded solid particles and planet formation.
Methods. We performed two-dimensional global simulations of the dynamics of gas and solids in a non-magnetized thin protoplanetary disk with the Pencil code. We used multiple particle species of radius 1, 10, 30, and 100 cm. We computed the particles' gravitational interaction by a particle-mesh method, translating the particles' number density into surface density and computing the corresponding self-gravitational potential via fast Fourier transforms. The dead zone is modeled as a region of low viscosity. Adiabatic and locally isothermal equations of state are used.
Results. The Rossby wave instability is triggered under a variety of conditions, thus making vortex formation a robust process. Inside the vortices, fast accumulation of solids occurs and the particles collapse into objects of planetary mass on timescales as short as five orbits. Because the drag force is size-dependent, aerodynamical sorting ensues within the vortical motion, and the first bound structures formed are composed primarily of similarly-sized particles. In addition to erosion due to ram pressure, we identify gas tides from the massive vortices as a disrupting agent of formed protoplanetary embryos. We find evidence that the backreaction of the drag force from the particles onto the gas modifies the evolution of the Rossby wave instability, with vortices being launched only at later times if this term is excluded from the momentum equation. Even though the gas is not initially gravitationally unstable, the vortices can grow to
in locally isothermal runs, which halts the inverse cascade of energy towards smaller wavenumbers. As a result, vortices in models without self-gravity tend to rapidly merge towards a m=2 or m =1 mode, while models with self-gravity retain dominant higher order modes (m=4 or m=3) for longer times. Non-selfgravitating disks thus show fewer and stronger vortices. We also estimate the collisional velocity history of the particles that compose the most massive embryo by the end of the simulation, finding that the vast majority of them never experienced a collision with another particle at speeds faster than 1 m s-1. This result lends further support to previous studies showing that vortices provide a favorable environment for planet formation.
Key words: accretion, accretion disks - hydrodynamics - instabilities - stars: planetary systems: formation - methods: numerical - turbulence
1 Introduction
The ill fate of the building blocks of planets in gaseous disks around young stars stands as one of the major unsolved problems in the theory of planet formation. Beginning with micron-sized interstellar dust grains, coagulation models predict growth to centimeter and meter size (Weidenschilling 1980; Dominik et al. 2007) in the denser environments of a circumstellar disk. Such bodies, however, are large enough to have already decoupled slightly from the sub-Keplerian gas, yet still small enough to be subject to a significant gas drag. The resulting headwind drains their angular momentum, leading them into spiral trajectories towards the star, on timescales as short as a hundred years at 1AU (Weidenschilling 1977a). Another acute problem is that such bodies have poor sticking properties and a low threshold velocity for fragmentation (Chokshi et al. 1993), such that collisions between them usually lead to destruction rather than growth (Benz 2000; Sirono 2004; Ormel & Cuzzi 2007). Such problems severely hinder growth to km-size by coagulation (Brauer et al. 2008a).
In view of these problems, other routes for breaching the meter size barrier have been pursued. A distinct alternative is gravitational instability of the layer of solids (Safronov 1969; Lyttleton 1972; Goldreich & Ward 1973; Youdin & Shu 2002). When the dust aggregates had grown to centimeter and meter size the gas drag is reduced and the solids are pushed to the midplane of the disk due to the stellar gravity. Although such bodies do not have enough mass to attract each other individually, sedimentation increases the solids-to-gas ratio by orders of magnitude when compared to the interstellar value of 10-2. It was then hypothesized (Safronov 1969) that due to the high densities of this midplane layer, the solids could collectively achieve critical number density and undergo direct gravitational collapse. Such a scenario has the advantage of occurring on very rapid timescales, thus avoiding the radial drift barrier.
This picture is nonetheless simplistic, in the view that even low levels of turbulence in the disk preclude the midplane layer of solids from achieving densities high enough to trigger the gravitational instability (Weidenschilling & Cuzzi 1993). Even in the absence of self-sustained turbulence such as the one generated by the magneto-rotational instability (MRI; Balbus & Hawley 1991; Balbus & Hawley 1998), the solids themselves can generate turbulence due to the backreaction of the drag force onto the gas. Such turbulence can be brought about by Kelvin-Helmholtz instabilities due to the vertical shear present in the sedimented layer of solids (Weidenschilling 1980; Weidenschilling & Cuzzi 1993; Sekiya 1998; Johansen et al. 2006), or by streaming instabilities induced by the radial migration of solids particles (Youdin & Goodman 2005; Johansen et al. 2006; Paardekooper 2006; Youdin & Johansen 2007; Johansen & Youdin 2007). In the turbulent motion, the solids are stirred up by the gas, forming a vertically extended layer where the stellar gravity is balanced by turbulent diffusion (Dubrulle et al. 1995; Garaud & Lin 2004).
But if turbulence precludes direct gravitational collapse through sedimentation, it was also shown that it allows for it in an indirect way. As solid particles concentrate in high pressure regions (Haghighipour & Boss 2003), the solids-to-gas ratio can be enhanced in the transient turbulent gas pressure maxima, potentially reaching values high enough to achieve gravitational collapse. Numerical calculations by Johansen et al. (2007) show that this is indeed the case, with the particles trapped in the pressure maxima generated by the MRI collapsing into dwarf planets when the gravitational interaction between particles is considered. They also show that the MRI is not necessarily needed, since the weak turbulence brought about by the streaming instability itself can lead to enough clumping under certain conditions. Another way of achieving high enough densities for gravitational collapse of the solid layer was shown by Rice et al. (2004) and Rice et al. (2006), where meter-sized solids concentrate prodigiously in the spiral arms formed in marginally gravitationally unstable circumstellar disks.
Such models, however, ignored the possibility of fragmentation of particles upon collisions. As the turbulence enhances the velocity dispersion of solids, destructive collisions become more likely. Moreover, upon destruction, the smaller fragments are tightly coupled to the gas and therefore dragged away from the midplane (Johansen et al. 2008), reducing the effective amount of solid material available for collapse. Such problem is particularly severe in the high mass disks investigated by Rice et al. (2004) and Rice et al. (2006), where the typical speeds of the boulders upon encounters are comparable to the sound speed.
The fragmentation problem could be avoided if the accumulation of solids happened, for instance, within a protective environment where the collisional speeds are brought down to gentler values. Anticyclonic vortices (Marcus 1990) have been shown to favor planet formation (Barge & Sommeria 1995; Tanga et al. 1996; Bracco et al. 1999; Chavanis 2000; Johansen et al. 2004) since, by rotating clockwise in the global counterclockwise Keplerian flow, they enhance the local shear and induce a net force on solid particles towards their centers. Klahr & Bodenheimer (2006) further argue that anticyclonic vortices would be less turbulent than the ambient gas, which in turn would lead to velocity dispersions that are low enough to prevent fragmentation. Vortices in disks can be the result of the baroclinic instability (Klahr & Bodenheimer 2003; Klahr 2004; Petersen et al. 2007), the Rossby wave instability (Lovelace et al. 1999; Li et al. 2000; Li et al. 2001) or, perhaps, the MRI (Fromang & Nelson 2005).
In this paper, we focus on vortices generated by the Rossby wave instability (RWI), which is a global instability where azimuthal modes experience growth in the presence of local extrema of a quantity interpreted as a combination of entropy and potential vorticity. In the linear phase, the instability launches inertial-acoustic waves. The non-linear saturation is achieved when the Rossby waves break and coalesce into anticyclonic vortices. It was shown by Varnière & Tagger (2006) that a favorable profile of the entropy-modified vorticity naturally arises if the disk has a slow-accretion zone, such as in the layered accretion model of Gammie (1996). In this model, ionization is provided by collisions in the hot inner regions, and by cosmic rays in the outer disk where the column densities are low (a standard value for the penetration depth of cosmic rays is a gas column density of 100 g cm-2). Throughout most of the midplane, however, the temperatures are too cold and the column densities are too thick for ionization to occur either way. The result is that, when threaded by a weak magnetic field, the disk displays MRI-active regions in the ionized layers, and a MRI-dead zone in the neutral parts around the midplane (Gammie 1996; Miller & Stone 2000; Oishi et al. 2007). Matter flows towards the star due to the high turbulent viscosity of the MRI-active layers, but upon hitting the border of the dead zone, it reaches a region of slow accretion and the flow stalls. However, as the flow proceeds unabridgedly from the outer active regions, a surface density maximum forms, which triggers the growth of the RWI.
The implications of this scenario for planet formation were first
explored by Inaba & Barge (2006), who use the RWI-unstable dead zone model
of Varnière & Tagger (2006) to study the accumulation of solids inside the
Rossby vortices. They confirm that the vortices are efficient particle traps,
since the solids-to-gas ratio was raised by at least one order of magnitude,
modeling the solid phase of the disk as a fluid. Such
approximation requires that the size of the solid particles be much smaller
than the gas mean free path. Since in the Minimum Mass Solar Nebula (MMSN;
Weidenschilling 1977b) at 5.2 AU the particles subject to maximum drift have a size comparable to
the mean free path, the sizes that a fluid approach can handle correspond to too
strong friction, thus ultimately underestimating the trapping performance of the vortical motion. In Lyra et al. (2008b, hereafter LJKP08), we
took the works of Varnière & Tagger (2006) and Inaba & Barge (2006) one
step further by including gravitationally interacting centimeter and meter
size solids treated as Lagrangian particles. In that Letter, we showed
that the solids concentrated in the vortices triggered by the RWI rapidly
reach critical densities and undergo collapse into rocky planets. The
resulting burst lead to the formation of 20 rocky protoplanetary embryos
in the mass range 0.1-0.6
,
along with hundreds of smaller bodies
following a mass spectrum of power law
.
In this paper we further detail the method used in LJPK08, also presenting a number of new results. In the following section we present the dynamical equations, followed by an in-depth analysis of the vortices in Sect. 3. In Sect. 4 we analyze the formation and evolution of the protoplanetary embryos, focusing on stability against erosion (Paraskov et al. 2006; Cuzzi et al. 2008) and tides from the gas, which we identify as an important disrupting agent. In Sect. 6 we investigate the response of the RWI to effects not considered in the original analysis of Lovelace et al. (1999) and Li et al. (2000), and in Sect. 7 we present a discussion of the limitations of the model. A summary and conclusions are presented in Sect. 8.
2 The model
2.1 Dynamical equations
We work in the thin disk approximation, using the vertically integrated equations of hydrodynamics
In the above equations G is the gravitational constant,











The function
is the drag force by which gas and solids interact.
In Eq. (8),
is the internal density of a
solid particle,
its radius, and
its
velocity relative to the gas. CD is a dimensionless coefficient that
defines the strength of the drag force. We use the formula of Woitke &
Helling (2003) that interpolates between Epstein and Stokes drag
where


where



2.2 Initial and boundary conditions
In this paper, we use the Pencil Code
in Cartesian and cylindrical geometry. The cylindrical runs do not include the gravity of
the particles and were therefore
only used for tests or runs without particles, as will become clear in the next sections.
A Cartesian box was used for the production runs. The Cartesian box ranges
.
The resolution is
,
unless stated otherwise. The
cylindrical grid ranges
,
with 2
coverage in azimuth.
The density profile follows the power law
and
the sound speed is also set as a power law
.
These slopes
are chosen to cancel a radial dependency on the mass accretion rate
if the viscosity follows the Shakura-Sunyaev
parametrization,
.
The velocity field
is set by the condition of centrifugal equilibrium
We use units such that










The dead zone is modeled as static viscosity jumps following arc-tangent profiles
where r1=0.6 and r2=1.2 are the locations of the jumps and




For the solids, we use 105 or
Lagrangian numerical particles.
For a gas mass of
and the interstellar solids-to-gas
ratio of 10-2, each numerical particle therefore is a super-particle
containing (in the lower resolution case)
of material. We use particles of radii
,
10, 30, and 100 cm, as also used in
Lyra et al. (2009). For our nebula parameters, maximum drift occurs for particles
of 30 cm, as detailed in that paper.
The particles are initialized as to yield a surface density following the same power law as the gas density, and their velocities are initialized to the Keplerian value. We use reflective boundaries for the cylindrical grid and frozen boundaries for the Cartesian grid (Lyra et al. 2008a). Both use the buffer zone described in de Val-Borro et al. (2006) to damp waves before they reach the boundary. Particles are removed from the simulation if they cross the inner boundary.
3 Vortices
The trapping mechanism of vortices is not only due to
its being a high-pressure region, but mainly due to the
vorticity of the flow. In an anticyclonic vortex (cyclonic
vortices are destroyed by the Keplerian shear), the motion occurs in the
same sense as the local shear, i.e., the gas rotates clockwise. Therefore,
at the antistellar point the angular momentum is decreased with respect
to a non-vortical flow; and conversely increased at the substellar point.
As a result, the gas at the antistellar point is accelerated
inwards, while the gas at the substellar point is accelerated
outwards. A net centripetal force
towards the eye ensues. The streamlines of vortices (or vortex
lines) are a set of Keplerian ellipses with the same semimajor
axis but different eccentricities, being circular in the center and
more eccentric outwards (Barge & Sommeria 1995). We show contours
of
on the surroundings of one of the giant vortices, in the
upper left panel of Fig. 1. As the gas drags the particles, the particles
also revolve around the vortex eye. But because the gas-solids
coupling is not perfect, the particles lose angular momentum and sink
deeply towards the center. In the next subsections we describe some of the
properties of the vortices present in our simulations
![]() |
Figure 1:
Enlargement around one of the vortices in a snapshot at 20
orbits. Contours of |
Open with DEXTER |
3.1 Launching mechanism - the RWI
The vortices in LJKP08 are triggered by the Rossby wave instability (RWI), a case of purely hydrodynamical instability in accretion disks. Considering azimuthal perturbations to the inviscid Euler equations, Lovelace et al. (1999) and Li et al. (2000) find that instabilities exist when the following quantity has a local extremum
![]() |
(14) |
The quantity

![]() |
(15) |
where
![]() |
(16) |
is the epicyclic frequency and
![]() |
![]() |
![]() |
(17) |
![]() |
![]() |
![]() |
(18) |
are the radial length scale of the entropy and density variations, respectively.






![]() |
= | ![]() |
(19) |
![]() |
= | ![]() |
(20) |
= | ![]() |
which in turn led Lovelace et al. (1999) to interpret


3.2 How sharp need the jump be?
Although formally the instability is triggered by a minimum or
maximum of
,
ones finds in practice that the amplitude and
radial width of the pressure bump present critical values beyond which
instability does not occur. Li et al. (2000) find that typically,
a pressure variation of 10-20% over a length similar
to the scale height of the disk is sufficient to trigger the instability.
The threshold, however, is problem dependent, depending - among other things -
on the geometry of the pressure variation (step jump or Gaussian,
for example; see Li et al. 2000, for details). Due to this, we
use a more general criterion to assess the threshold of instability.
Li et al. (2000) note that the threshold of instability for the RWI is
always reached before the Solberg-Høiland criterion for stability of
axis-symmetric disturbances
is violated. Eq. (21), therefore, provides a conservative estimate of whether or not the RWI is excited. The Solberg-Høiland criterion is easily understood. In a pressureless disk, the condition

![]() |
(22) |
is the square of the Brunt-Väisälä frequency, associated with the oscillations of buoyant structures in the presence of an entropy gradient. Physically, Eq. (21) means that a mode that is unstable/stable to shear can be stabilized/destabilized by pressure gradients and vice-versa.
We measure the epicyclic and Brunt-Väisälä frequencies in a series of 1D simulations where we varied the width
of the viscosity jump
(Eq. (13)). We find that for locally isothermal simulations, the
Solberg-Høiland criterion is violated at the outer edge
for
,
which is slightly less than one scale
height. In these isothermal simulations, the criterion depends
almost solely on the epicyclic frequency because, as the temperature
does not rise with compression, the pressure does not change enough for N2 to go appreciably negative.
To assess the effect of non-barotropic behavior, we replace the locally isothermal flow by an isentropic one with an adiabatic equation of state
which means that we allow heating and cooling by compression and rarefaction only, excluding viscous heating and radiative cooling. In Eq. (23),



![]() |
Figure 2: Evolution of the pressure bump with different widths of the viscosity jump ( upper most panels; the dashed line is the Heaviside jump, for comparison). The violation of the Solberg-Høiland criterion ( lower most panels) is a conservative indication that the threshold of instability of the RWI was reached. As the width of the jump increases, the threshold takes increasingly longer time to be breached. For jumps smoother than two scale lengths, the threshold was not reached (up to 500 orbits). |
Open with DEXTER |
The results are illustrated in Fig. 2, where we plot the viscosity profile
(upper panel) for different widths ,
and the time evolution of the
density, pressure and temperature bumps. The lower panels measure if the
Solberg-Høiland criterion was violated. We find that the pressure bump
sharpens considerably compared to the isothermal case, due to the high temperatures associated
with the compression. The consequence is that the Solberg-Høiland criterion
is violated by viscosity jumps up to
,
i.e., 3 times broader than
in the isothermal simulations. In this non-isothermal case, it is mostly the
Brunt-Väisälä frequency that leads
to negative values. The
effect of increasing the width is mainly of slowing the evolution of the quantity
towards negative values. It only takes five orbits for
,
but it takes 350 orbits when
is increased to 0.1. In
Fig. 2 we state
in terms of the scale height H=0.05 r0. We
present a resolution study of vortex excitation in Appendix A. We also
address the issue of vortex survival in a non-static dead zone in Appendix B.
3.3 Steady-state dead zone
If no transport happens in the dead zone, matter can do little more than piling up there as the inflow proceeds from the active layers. However, the accumulation of matter cannot proceed indefinitely since, as matter piles up, the conditions for gravitational instability would eventually be met (Armitage et al. 2001). The gravitational turbulence that ensues (Lodato 2008) would therefore empty the dead zone as the excess matter accretes, thus re-starting the cycle.
However, local simulations show that the dead zone has some level of
residual turbulence. This happens because the turbulence on the active
layers induce small levels of Reynolds stress in the dead zone
(Fleming & Stone 2003). If the inertia of the midplane layer is
not too high (Oishi et al. 2007), this forced turbulence
can lead to moderate
values with non-negligible transport
.
Terquem (2008) shows that steady state solutions in 1D models exist in this
case, as the dead zone gets denser and hotter to match the condition of constancy of the mass accretion rate with radius,
.
In this case, the steady state will have an
viscosity value in the active layers and a lower
in the dead zone.
Vortex formation by the RWI requires the presence of a pressure maximum.
In our model, and that of Varniére & Tagger (2006) and Inaba &
Barge (2006), the pressure maximum comes about by stalling the
accretion flow in the border of the dead zone. There is no requirement
that the dead zone should have zero viscosity, just a viscosity significantly lower than that of the active regions. We tested different values of /
,
and found that changing it from 0 to 0.1 has little effect on the instability. For higher values, the Solberg-Høiland criterion takes
increasingly longer to be violated. For
,
the Solberg-Høiland criterion is violated after 60 orbits. We also notice that the steady-state dead zones of Terquem (2008; see Fig. 3 of that paper) have the surface density varying by more a factor of
10 over a few scale lengths at the inner edge. Such profiles violate the Solberg-Høiland criterion, so the RWI is expected to be excited in
those conditions as well.
4 Embryos
4.1 Drag force cooling and compactness
![]() |
Figure 3: 1D simulation of collapse of 1000 particles. With drag force the kinetic energy of the particles is efficiently dissipated and the particles collapse at subgrid scale towards infinite density. When the drag force is excluded the system cannot dissipate energy and a halo of particles, 10 grid cells wide, is formed. |
Open with DEXTER |
The embryos formed in our simulations present a number of interesting features.
We first would like to address the issue regarding their physical size. The
embryos consist of a cluster of a large number of particles, held together by
their collective gravitational pull. But are they strongly bound like solid
objects? Or do they consist of loosely coupled objects in the
same potential well? To answer this question we measure the rms spatial
dispersion of the particles inside the cluster, defined as
where n is the number of particles within the Hill sphere of the clump,



Such compactness is due to the dynamical cooling provided by the drag force. We illustrate this in Fig. 3. The figure shows the results of 1D simulation with a thousand interacting particles with and without gas drag. Without gas drag the particles have no means to dissipate energy and perform oscillations about the center of mass. The very inner particles show virialization, while the outer particles form a halo extending for a radius of 10 grid cells in average.
When including gas drag, the system gets so dissipative that the kinetic energy is soon lost and the ensemble of particles collapses. The exponential decay of the particles' rms position seen in the upper left panel of Fig. 3 shows no sign of flattening, down to a millionth of a resolution element. This leads us to infer that collapse to zero volume is ongoing. This is of course expected, since no mechanism to provide support against the gravitational pull is present.
In view of this, the question is why our planets, which are
subject to drag forces, do not shrink to zero as well but stabilize at a
very small but finite radius. We are drawn to the conclusion that this is
a numerical issue. The tests of Fig. 3 were done with a fixed
time-step. But when the particles cluster together to form a planet in
our simulations, they end up dominating the time-step. The position
update
therefore occurs
with the maximum
allowed by the Courant
condition, which is that the fastest particle should move by one grid cell.
Due to this, the time resolution of the subgrid motion around the center of
mass of the cluster is under-resolved. With the overshot
,
the particles that are attracted towards the center
of mass of the clump will end up in a position past it. In the next time
step they will be attracted to the center of mass from the other side, but will
once again overshoot it. The result of this is that the particles will execute
undamped oscillations, leading to a finite rms radius. We performed
tests like those of Fig. 3 with a particle-controlled variable
timestep, confirming this explanation. We conclude that the fact that the most
massive embryo has a stable rms radius compatible with its mass is but a deceptive
coincidence.
We stress that this drag force cooling will cease to be efficient as
the solids-to-gas ratio grows too large (
), because in
this case the backreaction would be too strong and the gas would simply be
dragged along with the particles. In this case, a Keplerian disk of solids
might form, accreting matter onto the planet due to collisions between the
orbiting solids. This accretion regime is nevertheless beyond the scope of the
current paper.
4.2 Mass loss
In Fig. 2a of LJKP08 we showed the evolution of the most massive clump.
However, as the clumps have different feeding rate and some of them experience
mass loss, the most massive clump shown there is not always the same clump.
In Fig. 4 we contrast this with the evolution of the
most massive clump at the end of the simulation, which we tracked backwards
in time. Such clump started in the inner disk, showing 0.8
by 40 orbits.
By this time, the most massive clump was a 3
clump in the outer disk.
![]() |
Figure 4:
Evolution of the most massive clump formed (solid line), which we traced back in time from the end of the simulation. It differs from the instantaneous most massive clump (dashed line) because the clumps have different feeding rate and can also experience mass loss, as in the episode that happened at |
Open with DEXTER |
The most remarkable feature of this plot is the mass loss event at 90 orbits.
Figure 5 shows that it consists of the detachment of a 0.8
object from the original cluster, of 6.7
.
The detachment is already
seen at 87 orbits, although the separation is quite small (4 times the Earth-Moon mean
separation,
). At 89 orbits, the smaller object left the Hill sphere of the more massive embryo. They finally appeared as separate objects, and the maximum mass decreased.
![]() |
Figure 5: Time series of the mass loss episode. Due to gravitational tides from the massive vortex whence the embryo formed, a large chunk of particles was detached from the original cluster. At 87 orbits, a separation 4 times the Earth-Moon distance is seen. The separation grows and two orbits later the two bodies do not overlap Hill radii, thus counting as separate objects. |
Open with DEXTER |
We see evidence that this puzzling behavior is due to gravitational tides from the gas. The gas is too pressure-supported to undergo collapse, but the vortices concentrate enough material to yield a considerable gravitational pull. We illustrate this in Fig. 6, where we show the state of the disk before the mass loss episode (at 80 orbits, Figs. 6a-d) and after that (at 100 orbits, Figs. 6e-h). The plots show the surface densities of gas and solids, and the potential associated with them. Even though the clumping of solids yield a considerable gravitational pull (Figs. 6d, h), these figures show that the dominant contribution to the gravitational potential comes from the gas - more specifically from the vortices, where the gas density peaks one order of magnitude denser than the initial condition.
The most massive clump is located in the inner disk at (X,Y)=(-0.40, -0.53) in Fig. 6b, not clearly identifiable amidst the other particles trapped inside the vortex. However, the embryo is immediately observable as the bright point at (X,Y)=(-0.65, -0.19) in Fig. 6h (also visible in Fig. 6f, albeit less prominently). There are two features in this plot that are worth noting. First, by comparing the location of the embryo in these figures with the location of the vortices, we notice that the planet has left its parental vortex. Second, the inner vortices have undergone the transition from the m=3 to the m=2 mode. Due to merging, their gas density has increased, with dramatic consequences for the embryos within them.
![]() |
Figure 6: The state of the disk before a)- d) and after e)- h) the mass loss episode. The conspicuous difference between them is that the embryo has left its parental vortex from one snapshot to the other. It is seen as a bright spot in panels f) and h), at (X,Y)=(-0.65, -0.19). In panel b) (before the mass loss), the embryo is at (X,Y)=(-0.40, -0.53) but not easily spotted among the swarm of solids inside the vortex. Panel i) shows a horizontal slice through this location, in which we see that the density of solids does not peak much higher than the gas density at the location of the embryo (panel j)). Significant gas tides are expected as the gravitational potential (panel k)) and acceleration (panel l)) have similar contributions from the gas and solid components. |
Open with DEXTER |
We assess how the gravity of the gas influences the motion of the particles (Fig. 6i-l). In Fig. 6i we take a horizontal density slice at the position of the most massive embryo at 80 orbits. Fig. 6j is an enlargement of Fig. 6i around X=-0.53, where the embryo is located. We see that the densities of solids and gas peak at similar values. The next figures show the gravitational potential (Fig. 6k) and acceleration (Fig. 6l) around the embryo. The gas produces a deeper gravitational well, albeit smoother than the one displayed by the solids. In the acceleration plot it is seen that the pull of the gas is greater than the pull of the embryo already at a distance of just 0.26 AU (0.03 in code units, corresponding to two grid cells) away from the center. And even where the pull of the solids is strongest (one grid cell away from the center), the gravity of the gas still is an appreciable fraction of the gravity of the solids. Tides from the gas are unavoidable.
It is beyond the scope of this paper to consider the full
mathematical details of the theory of tides, especially because the two bodies (the vortex and the
embryo) are extended. Instead, we consider the following toy model. The tidal force
experienced by the planet is proportional to the gradient of the acceleration ainduced by the gas. It is also proportional to the radius R of the planet:
.
Since
,
according to the Poisson equation, the tidal force is proportional to the local value of the density
![]() |
(26) |
We consider the 3D volume density to avoid the requirement of using the Dirac delta in the 2D case. Considering the planet spherical, Newton's second theorem holds and we can write FG = -GM/R2 for the planet's (self-)gravitational force at its surface. Substituting


i.e., the ratio of the disrupting tidal stresses to the self-gravitating forces that attempt to keep the planet together is directly proportional to the gas-to-solids ratio. At 80 orbits, as seen in Fig. 6j this ratio is around unity. As the vortices undergo the transition from the m=3 to the m=2 mode, their peak density increases (while the planet remains at a constant mass), the tides eventually become strong enough to cause the mass loss event of Fig. 5. This effect will probably be less dramatic in 3D simulations because, as the particles settle in the midplane, the ratio of the volume gas density to the bulk density of solids



4.3 Erosion?
Cuzzi et al. (2008) points that erosion is of prominent importance in the
stability of self-gravitating clumps of particles. They put forth a model where self-gravity plays a role analogous to that of surface tension in liquid drops, preventing disruptionagainst ram pressure forces from the gas. The clumps are held together by self-gravity if the gravitational Weber number (in analogy with the surface tension case) is less than a critical value, close to 1. The gravitational Weber number is defined as the ratio of the drag to self-gravitational accelerations
Cuzzi et al. (2008) further point that in numerical models, artificial viscosity can largely exaggerate the disrupting effect of the ram pressure. This happens because, as the clumps are small, they are deeply within the viscous range of the grid, whereas in the real solar nebula the dissipation happens at much smaller scales. The Reynolds number of the flow past the particles is therefore much smaller than what a real clump would experience
![[*]](/icons/foot_motif.png)
It is interesting to assess if this mechanism plays a significant role in our models, or even if it can account for at least some of the mass loss events.
We can estimate how important erosion will be for our clumps the following way.
We approximate the clumps as single point masses so that
,
where M is the total mass of the clump. Plugging this in Eq. (28), we write the gravitational Weber number as
![]() |
(29) |
For Epstein drag (Eq. (10)), CD does not depend on



We can simplify the
by writing the mass M as
,
and the drag force as
.
Thus, at
,
We confirm in Appendix C that Eq. (30) is sufficiently accurate in predicting the onset of erosion. For our choice of parameters,
so for a flow of







4.4 Can the embryos be held together indefinitely?
To answer the question of how long lived these clumps are, we take the last snapshot of the simulation and switch off the hydrodynamics. The particles now move under the influence of the stellar gravity and their own self-gravity only. We run for additional 200 orbits to assess is self-gravity alone can maintain the cluster of particles together. The result is shown in Fig. 7.
The clustered particles do not disperse, and the most massive embryo maintains the same mass for the additional simulating time. We do not see difference in the degree of compactness of the most massive embryo. As the vortices are shut down, the unbound 1 cm sized particles that were too well coupled to the gas to be dragged into the eye are released from they vortical confinement and spread over a wider annulus.
We have no reason to suspect that the situation will change over longer timescales. We conclude that the embryos can be held together indefinitely.
4.5 Collisions
As stated before, one of the problems that solids accumulation inside vortices can potentially solve is the issue of fragmentation of particles upon collisions, a drawback for both coagulation (Brauer et al. 2008a) and gravitational instability models (Rice et al. 2006; Johansen et al. 2007) of planetesimal growth.
To assess if fragmentation poses a significant barrier for the formation of the protoplanetary embryos in this study, we take the most massive embryo by the end of the simulation, flag the particles that constitute it and trace them back in time, calculating their collisional velocity history. The collisional speed for each particle is calculated by taking the closest neighbor to that particle within the range of a grid cell. There is, however, a subtlety concerning the difference between collisional velocities and relative velocities. A collision between particles i and j only happens when the separation rij between them tends to zero. For our resolution and choice of r0, a grid cell is 0.08 AU wide, thus existing plenty of room for subgrid Lagrangian dynamics. In particular, the velocity difference due to the Keplerian shear between the inner and outer radial borders of a grid cell can be significant. At the inner edge of the dead zone of the model presented in this paper (3.12 AU) for instance, this difference amounts to 434 m s-1. As this velocity difference is due to the separation between particles, it vanishes when rij tends to zero, thus never contributing to the true collisional speed. In Fig. 8 we show these uncorrected relative velocities in the initial condition and in a snapshot at 40 orbits, plotted against separation. The clear correlation follows what is expected from the Keplerian shear (dashed line) in the initial condition.
![]() |
Figure 7:
Can self-gravity alone maintain the formed embryos together? We switch off the hydrodynamics at the last snapshot and run for additional 200 orbits. The most massive embryo, which had already left the parental vortex (at |
Open with DEXTER |
![]() |
Figure 8:
The relative velocities between particles are contrasted with
collisional velocities calculated from Eq. (32). At the
initial condition (t=0) the relative velocities between particles follow
what is expected from the Keplerian shear (dashed line), reaching as far
as 300 m s-1 inside a grid cell (which is 31
|
Open with DEXTER |
The gas motion adds another velocity that has to be taken into
account. Solid particles are dragged by the gas motion, yet gas
streamlines never intersect. The gas motion itself
thus introduce velocities that never participate in collisions. We
correct for these two by the following procedure. For each particle
in the pair involved in a collision, we consider its velocity
relative to the gas (the same quantity that appears in the
drag force,
). We then define the
collisional velocity vector as
Equation (32) ensures that tracer particles, which follow the gas streamlines, should never experience collision. Furthermore, as the shear dependence of



The results of the collisional velocity history of the particles that constitute the most massive embryo at the end of simulation are plotted in Fig. 9.
Figure 9a shows a cumulative plot of the mass
of the embryo, which defines the time t0 that each flagged particle
was accreted. We define the time of accretion as the moment when the
particle approached the grid point
nearest to the maximum
of particle number density (also defined by the flagged particles) by
less than
(the grid cell diagonal) and kept
until the end of the simulation. Although this is not as strict as the definition of accretion we have been using before (based on the Hill criterion and escape velocity), this simpler criterion captures what happens before collapse (i.e., before the maximum of particle number density becomes a bound protoplanetary embryo) and serves well our purpose of illustrating the behavior of collisional speeds at close separations. The first episode of accretion takes place at 40 orbits, coinciding with the time that the clump of particles became bound (in accordance with Fig. 4). Further accretion proceeds over the next 18 orbits, with the maximum mass being attained at 58 orbits. No other particle was accreted after this time.
Figure 9b shows
histograms of the maximum collisional speed that a particle experienced
before (dot-dashed line) and after (solid line) accretion. The latter
refers to the maximum of
taken between t0
and the end of the simulation (200 orbits). The former refers to a
time interval of 10 orbits before t0. As the vortices in the inner
disk just fully develop at
30 orbits, the dot-dashed histogram
is more representative of a situation where particles are inserted in a
disk with existent vortices. The conclusion is striking: the vast
majority of the particles that constitute the embryo never experienced
a collision more violent than 1 m s-1.
![]() |
Figure 9:
a) Accretion history of the most massive embryo.
Collapse happens at 40 orbits, and further accretion of particles happens through the next 18 orbits, after which the maximum mass is attained. This plot defines the time t0 that each particle is accreted. b) Histograms of the maximum collisional velocities experienced by a given particle, before (dot-dashed line) and after (solid line) its accretion time t0. The solid line histogram represents the maximum collisional speed from t0 until the end of the simulation. The dot-dashed line histogram covers a time interval of ten orbits before t0. As the vortices form at |
Open with DEXTER |
In Fig. 9c we show the maximum collisional velocities
of Fig. 9b as a function of the separation between a given
particle and its closest neighbor, also before and after accretion.
The distribution before accretion is trimodal, with particles with low speeds (<1 m s-1) at small separations, particles with high speeds (<20 m s-1) at small separations, and particles with high speeds at large separations (>10
). Only the second group of particles would have undergone fragmentation. The first group is below the fragmentation velocity threshold, whereas the large separations of the third group imply they never experienced an encounter close enough to lead to a collision. After accretion, virtually all particles (99%) belong to the first group.
In Fig. 10 we plot the collisional velocities versus separation at selected snapshots instead of historical maxima. In these plots, we only used the particles for which Eq. (33) was satisfied, i.e., considering only the collisions that are participating on the formation process of the embryo. At 30 orbits, a small number of particles is observed (87), 78% of these showing safe collisional speeds (<10 m s-1). 5 orbits later 119 particles are within the grid cell of the forming embryo, 92% of which show gentle collisions. At the time the overdensity gets bound (40 orbits), it is formed by 639 particles, with just 7 of these showing collisional speeds greater than 10 m s-1. At 45 orbits, all 877 particles display safe speeds. The tendency seen in this time series towards an increasing number of encounters at low separations and at low collisional speeds indicates that collapse towards zero volume is ongoing. Indeed, at 70 orbits, we observe that most of the particles occupy the same point in space.
5 Size distribution
![]() |
Figure 10:
Collisional speeds as a function of separation at selected snapshots at 30, 35, 40, and 45 orbits (
|
Open with DEXTER |
In our simulations, we considered the solid phase of the disk represented by particles of 1, 10, 30, and 100 cm. In this section we discuss a number of issues, relevant to the simulations, related to having a size spectrum instead of single-phasing.
5.1 Aerodynamical sorting
One of the most prominent features of the embryos formed in our models is
that they are composed primarily of same-sized particles. This is mostly
due to aerodynamical sorting. As particles of different size have different
friction times, differential drag occurs inside the vortex, effectively
sorting the particles spatially by size. Moreover, the stationary point
is determined by a balance between the Coriolis and the drag force, in
such a way that the eye of the vortex is the stationary point only
for
,
or perfect coupling. In general, the stationary point is azimuthally shifted with respect to the eye, according to the particle size (Youdin 2008).
The aerodynamical sorting inside the vortex can be seen in Fig. 11, which corresponds to the vortex of Fig. 1. As similar particles drift alike, streams of same-sized particles are clearly seen in the vortical flow.
![]() |
Figure 11: Aerodynamical sorting for the particles trapped in the vortical motion. The figure is centered at the vortex shown in Fig. 1, at 18 orbits. The unit of length is the Earth-Moon mean distance, and the Y-coordinate points to the star. As particles of different size have different friction times, differential drag occurs inside the vortex, effectively sorting the particles spatially by size. |
Open with DEXTER |
5.2 Differences between the inner and outer embryos
![]() |
Figure 12:
The radius of the particle subject to maximum drift (
|
Open with DEXTER |
Another interesting feature of our results is that once the vortices are formed, they easily trap the 10 cm and 30 cm particles both in the inner and in the outer edge of the dead zone. Yet, by the end of the simulation there seems to be a preference for the embryos in the inner disk to be composed of larger particles (30-100 cm), while in the outer disk, more embryos formed of the smaller particles (1-10 cm) are seen. We explain these two features in the following sections.
5.2.1 Preferential sizes in different locations of the nebula
The first feature (inner and outer vortices equally trapping particles of same size) follows from the fact that the general drag coefficient (Eq. (9)) yields a nearly flat profile for the radius of the particle with maximum drift (
)
versus distance. This is seen by calculating the stopping times
of the different particles as a function of the dynamical variables
![]() |
= | ![]() |
|
= | ![]() |
(34) |
and calculating the radius for a given stopping time





![]() |
(35) |
and, as the radius is positive, the solution is unique
![]() |
(36) |
Figure 12 shows this curve as a function of radius for



5.2.2 Tidal disruption and erosion of the inner embryos of
10 cm
![]() |
Figure 13: Time series of the gas and solid density (for each individual particle size), showing the destruction of the embryos of 10 cm particles when the inner vortices undergo the transition from the m=3 to the m=2 mode. At 75 orbits ( upper panels), embryos are seen in both the inner and outer vortices, for both 10 cm and 30 cm particles. At 85 orbits one of the embryos in the 10 cm phase was disrupted, followed by a second at 95 orbits, and the last one ten orbits later. The embryos composed of 30 cm particles also experience tides. But as their gravitational Weber numbers are smaller, they just undergo splitting, the large fragments being more stable against erosion than the embryos composed of particles of 10 cm. |
Open with DEXTER |
![]() |
Figure 14: Upper panels: evolution of a disk without particles and without self-gravity, which serves as a control run for the next plots. Middle panels: evolution of a disk without particles but with self-gravity. The difference compared to the upper panels is that the dominant mode in the outer disk changed from m=2 to m=5. Self-gravity modifies the dispersion relation of the RWI, or it stalls the inverse cascade of power known to occur in 2D turbulence, or both. Lower panels: evolution of a disk without self-gravity but with particles. The backreaction leads to an early excitation of the RWI in the inner edge of the dead zone. Conversely, the outer edge goes unstable later when compared to the other two runs. Since the particle density is not high enough to excite the streaming instability, we take it as evidence that the backreaction modifies the dispersion relation of the RWI. |
Open with DEXTER |
We see evidence that the second feature (the absence of clumps of 10 cm in the inner disk at later times even though they were formed) is due to tidal disruption in the same episode that lead the most massive embryo to lose mass. We illustrate this in Fig. 13, a time series of the gas and the solid phases, the latter split into the 4 different particle sizes. The upper plots, at 75 orbits, show that embryos of 30 cm and 10 cm were formed in both the inner and outer vortices.
The gas plots of the time series of Fig. 13 illustrate the transition from the m=3 to the m=2 undergone by the inner vortices, as mentioned in Sect. 4.2. This raises their density so the tides get stronger. The plots of the 30 cm phase show that the embryos composed of these particles split into smaller objects, which nevertheless can still keep their physical integrity.
The fate of the clusters composed of particles of 10 cm is different, though. As
the gas density increases, so does the gravitational Weber numbers of the embryos.
Erosion starts to play a more significant role. As the density increases inside the
inner vortices (a factor 5 relative to the initial condition at 75 orbits; 8 at 200 orbits), the embryos of 10 cm particles start to behave more and more like the
clusters of Fig. C.1. At 85 orbits, one of embryos
of 10 cm particles in the inner edge was destroyed. At 95 orbits, a second
embryo was disrupted. At 105 orbits, the third embryo was also destroyed
by the combined effect of tides and erosion. 25 orbits later, the 10 cm particles
have dispersed through the inner edge of the dead zone. The tides from the gas prevent
them from assembling once again.
The outer vortices never get as strong as the inner ones. The result is that although the inner embryos of 10 cm particles are destroyed, the outer ones are kept until the end of the simulation.
6 The response of the RWI to the drag backreaction and self-gravity
The evolution of the RWI was studied analytically for the case of a low mass dustless disk only. In Fig. 14 we show how the effects of gas self-gravity and backreaction from the solids affect the evolution of the instability.
The upper panels of Fig. 14 show a disk without solids and without self-gravity. In the middle panels we included self-gravity, while in the lower panels we included solids. The appearance of the disk in the three simulations is shown in selected snapshots at 5, 10, and 15 orbits.
The self-gravitating and non-selfgravitating dustless cases (upper and middle panels) look similar, with the RWI being excited first in the outer edge of the dead zone. However, there is a crucial difference between them. The snapshot at 15 orbits shows a prominent m=2 mode in
the outer disk for the non-selfgravitating case, while the run with self-gravity displays a dominant m=5 mode at the same time. This puzzling result is made even more interesting by recalling that the gas is gravitationally stable (
). That the growth rates of different modes vary that significantly for such a value of Q is indicative that the dispersion relation of the RWI is probably remarkably sensitive to self-gravity.
The simulation with drag backreaction (lower panels) also displays a number of interesting features. First, the RWI was excited in the inner disk as early as 5 orbits. In contrast, the control run
without solids (upper panels) has the instability appearing first in the outer disk, and at later times (10 orbits). The conclusion is that the particles induce vorticity on their own. Even though it is clear that this behavior has to do with free energy being transfered from the particle motion to the gas motion, it is not obvious if this result can be linked to the streaming instability (Youdin & Goodman 2005) since the solids-to-gas ratio is not nearly as high as the one needed to excite it (
). Instead, it is more likely that the backreaction
is modifying the dispersion relation of the RWI.
Another interesting feature of this run is that although the RWI was excited in the inner edge of the dead zone as early as 5 orbits, the outer edge just went unstable as late as 15 orbits. In contrast, the control run (upper panels) shows the outer disk going unstable at 10 orbits. As the backreaction hastens the growth of the RWI in the inner edge, it is unclear why it should stall it in the outer edge. One possibility is that the Rossby waves launched by the edge that first goes unstable interferes destructively with the perturbations fighting to grow in the other edge.
The dominant mode also changed from the dusty to the dustless case. The latter has m=4 and m=5 modes being dominant in the inner edge. In the dusty case a m=2 mode is seen instead. However, since in the dustless case it is the outer edge that displays a m=2 mode, another explanation comes to mind. As the models are 2D, we are probably witnessing the inverse cascade phenomenon due to enstrophy conservation. The vortices are simply cascading energy towards the larger scales, so the edge that goes unstable first (outer in the dustless case, inner in the dusty) will also reach a dominant m=2 mode first (possibly also m=1 at later times).
![]() |
Figure 15: In the simulation shown in LJKP08, the vortices maintained a m=4 mode in the outer edge of the dead zone until the end of the simulation. By switching self-gravity off, we see that in less than 15 orbits the outer vortices merged into a m=3 mode, while the m=3 mode in the inner disk edge into a m=2 mode. This is because as Q decreases, the size of the vortices approaches the Jeans length scale, which effectively halts the inverse cascade of energy. |
Open with DEXTER |
If this is the case, then self-gravity somehow halts the inverse cascade that took place in the dustless non-selfgravitating case. It is also instructive to compare the dusty non-selfgravitating run (lower panel) with the dusty selfgravitating run of LJKP08 (Fig. 1 of that paper). In that case, the m=4 was dominant in the outer edge of the dead zone until the end of the simulation at 200 orbits. We also perform a test (Fig. 15) that consists of switching the self-gravity off in the run of LJKP08, and checking the evolution of the vortices. Without self-gravity, the m=4 mode turns into a m=3 mode in less than 15 orbits. In the inner edge, a m=2 mode developed out of the otherwise dominant m=3 mode. The inverse cascade indeed resumed.
This result was also very recently reported by Mamatsashvili & Rice (2009). Without self-gravity, the vortex size is limited by the pressure scale height H. Once vortices grow to sizes of a few times H, the vortical flow becomes super-sonic. The vortex then radiates density waves that carry energy away and limit further growth. Mamatsashvili & Rice (2009) point that in the presence of self-gravity, the Jeans length
,
where
is the
Toomre Q parameter, poses another limitation to the maximum
size of a vortex. In Fig. 15 we see that the vortices in
the self-gravitating runs show
.
The growth seen when
self-gravity is switched off is a result of this constrain being lifted.
7 Limitations of the model
The presented models are admittedly simplified. In this section, we state what we consider the main limitations of our calculations to be.
7.1 Two-dimensionality
The most stringent limitation of the models is the 2D approximation, which leads to a number of features, stated below.
7.1.1 Vortex formation and survival
The question of the excitation and sustainability of vortices in three
dimensions is the matter of an old, yet unsettled, debate. Once excited,
anticyclonic vortices are easily maintained in 2D simulations where, unless
viscosity is present, they cannot decay and will instead merge, growing in
size in a cascade of energy towards the largest scale of box (e.g., Johnson &
Gammie 2005). However, three-dimensional studies in the context of
protoplanetary disks found that tall vortex columns are destroyed, both
in non-stratified (Shen et al. 2006) and in stratified (Barranco & Marcus
2005) local boxes. This phenomenon is understood as a result of the elliptic instability (Crow 1970; Gledzel et al. 1975; Kerswell 2002), by which the stretching term
,
absent in 2D, breaks down elliptical streamlines such as vortical flow. For a vortex to grow in 3D, the baroclinic term
has to counter the stretching term.
An indication that vortices can be sustained in three dimensions is present in the study of Edgar & Quillen (2008). These authors simulate a stratified disk in spherical coordinates with an embedded giant planet. In their inviscid run, the RWI is excited, leading to Rossby vortices at the edges of the gap, much like as in the 2D runs of de Val-Borro et al. (2007). The vortices launched in three dimensions are long lived and vertically extended, apparently following the same scale height as the surrounding disk. We remark that the MRI-generated vortex of Fromang & Nelson (2005) is also seen to be long-lived in a unstratified global disk. The studies of Edgar & Quillen (2008) and Fromang & Nelson (2005) both use a locally isothermal equation of state, which has large-scale non-zero baroclinity due to the static radial temperature gradient. Furthermore, the existence of the RWI in 3D is demonstrated by the simulations of Méheut et al. (2008).
7.1.2 Strength of the vortices
A major impact of the 2D approximation is the inverse cascade due to enstrophy conservation that overpowers the vortices. An in depth study of the formation, development and structure of Rossby vortices in 3D global accretion disks is needed to realistically address the issue of planet formation inside these structures.
7.1.3 Particle sedimentation
Another limitation posed by the two-dimensionally is that the particles and the
gas have the same infinitely thin scale height. The result of this is that
the back-reaction of the drag force from the particles onto the gas is
largely underestimated in
our models. In 3D disks, the midplane particle layer is far denser
due to sedimentation, so the ratio
is far greater than the
ratio
used in Eq. (3). The stronger
backreaction that ensues is known to excite the streaming instabilities if
(Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007). This instability enhances particle clumping, thus aiding collapse (Johansen et al. 2007). However, the effect of this strong backreaction on the vortices is poorly known.
7.1.4 Different particle scale heights
As the particles sediment, what sets the particle scale height is the equilibrium between turbulent diffusion and vertical gravity. Controlled by the drag force, the turbulent diffusion depends on the particle radius, and so does the equilibrium scale height of the solids (Dubrulle et al. 1995). Because of this, particles of 1 m radius settle in a thinner layer than those of 1 cm particles. Inside a vortex, turbulent motions are expected to be weaker (Klahr & Bodenheimer 2006), bringing the layer of solids closer to a 2D configuration, but a dependence on radius is still expected. We could not model this effect on our simulations.
7.1.5 Gas tides and mass loss
The strength of the disrupting gas tides is yet another effect related
to the difference between 2D and 3D models. As discussed in Sect. 3.2, the
tides are proportional to the gas-to-solids ratio
,
thus expected to be much weaker in 3D where sedimentation considerably increases
in the midplane. As the vorticity is also expected to be weaker in 3D, the peak of
at the vortex's eye would be weaker than
in a 2D calculation, further weakening the effect of tides. We therefore anticipate
embryos formed in 3D simulations to be significantly less prone to mass loss
than the ones presented in this study.
7.2 Equation of state
In this study, we used very simple equations of state: isothermal
(Eq. (7)) or adiabatic (Eq. (24)). The effect of the equation
of state can be appreciated by seeing the evolution of the Solberg-Høiland criterion
in isothermal and adiabatic simulations. In the former, it is the epicyclic
frequency that brings
to negative values, while in the latter the
criterion is broken mainly by the Brunt-Väisälä frequency.
The excitation of the Rossby wave instability is greatly favored in the presence
of a strong entropy gradient, and made more difficult (yet not impossible) as the
disk approaches isothermality. Therefore, an accurate modeling of the energy
budget - solving for radiative cooling and turbulent heating -, is something to
pursue in order to more realistically address the evolution of the RWI and the issue
of planet formation inside the vortices that constitute its saturated state.
7.3 Aerodynamics of the embryo
The aerodynamics of a super-particle is controlled by the radius
of
the individual rocks. This means that even though the ensemble of
rocks has the same position and velocity, there is still space between them and therefore they have contact with the gaseous nebula through all their surface area.
However, after gravitational collapse occurs, the solids are not any longer an ensemble of pebbles and boulders with free space between them, but a single massive object of planetary dimensions. This leads to a radical change in the aerodynamical properties. Yet, in our simulations, we still consider the collapsed body as an ensemble of super-particles, with the aerodynamical properties of individual pebbles and boulders. This is certainly a limitation of the model.
For a large planet, the correct treatment would be to consider that, after collapse, we leave the regime of particle-gas Stokes drag and enter planet-disk interaction by gravitational friction (type I migration). In the solar nebula the two drags have similar strength for bodies of 100 km. As we solve for the self-gravity of the gas, the latter is included in our model, albeit limited by the resolution of the grid. The fact that we keep using Epstein-Stokes drag on the super-particles after collapse might make an embryo more stable, especially in view of the very effective dynamical cooling provided by the drag force (Fig. 3). Substituting collapsed clusters by a sink particle that feels the gas gravity but not the gas drag is a possible solution, but also has caveats on its own. The evolution of sink particles depends too much on artificial numerical parameters such as the accretion radius. Furthermore, a sink particle does not suffer tidal effects, which we showed to be non-negligible.
7.4 Coagulation and fragmentation
As dust particles are drawn together, electromagnetic interactions occur at their surface, causing sticking under favorable conditions. Brauer et al. (2008b)
show that density enhancements like the ones we see - where matter accumulates due to a discontinuity in viscosity -, dramatically
favor coagulation. As particles are drawn together and the relative velocities are
reduced, growth by coagulation occurs for a range of mass accretion rate
and
threshold fragmentation velocity
.
They find that the meter size barrier
can be breached for mass accretion rates up to
(for
m s-1) and threshold fragmentation velocities no
less than
m s-1(for
).
This raises the possibility that even before the RWI excites the vortices, coagulation will have depleted the population of centimeter and meter sized objects onto bodies that are too large for our proposed mechanism to occur efficiently. As we see, it is preferentially the 10 and 30 cm sized particles that concentrate into planetary embryos.
The timescale for coagulation, however, is much longer than the time-scale
for gravitational collapse. We see growth to Mars size taking place
in only five orbits (60 yr), while growth by coagulation
from meter to kilometer size occurs on timescales of a few thousand years
according to Brauer et al. (2008b). On the other hand, it could as
well be that the favorable environment provided by the vortices act as to
speed up coagulation even further. This, of course, is not bad. Growth beyond
the preferred size will lead to decoupling from the gas and ejection from the
vortex that, in the end, behaves as a planetesimal factory.
Fragmentation is an important piece of reality that we did not include in our
model. Nevertheless, we showed in Sect 4.5 that the majority of the particles
were never involved in collisions with speeds greater than 1 m s-1. These are
of course very good news for planet formation. However, we feel the need to
stress that the time interval between snapshots in Fig. 9a is of whole orbits (totaling 200 snapshots). The number of high-speed
impacts could be greater had we checked the collision speeds at every
time-step instead. Although desirable, this would have made the simulations
computationally very expensive since it must be done in runtime. The result
of Fig. 9b should therefore be taken only as further evidence
for low collisional speeds inside vortices, not as conclusive proof of it.
Carballido et al. (2008) further point that the low collisional
speeds at low separation may be unrealistic. This is because the particles
couple to the smallest eddies, whose size is a function of the mesh
Reynolds number. These authors find average collisional speeds of 0.05
for particles of stopping time
(
cm
in our models), but notice a sharp decrease of the collisional
speeds towards smaller separations. In our simulations, we are considering
encounters that happen inside a grid cell, where we do not resolve the
velocity field, so this may indeed be the reason behind the low collisional velocities
we find. However, we point that there is a difference between the simulations of
Carballido et al. (2008) and those presented in this paper. They considered
particle concentrations in the transient pressure maxima of the
turbulence generated by the MRI, whereas we consider particle concentrations
within long lived anticyclonic vortices. As vortex structures tend to
merge and grow, there is less power available at the smallest scales when
compared to MRI turbulence.
7.5 Disk mass
We stated in Sect 2.2 that the constrain of having 30 Earth masses of solid in the inner disk, together with the interstellar dust-to-gas ratio of
,
and a surface density slope of -0.5 implies 0.01 solar mass within the modeled range. This disk can be considered
massive since the shallow slope then implies a massive outer disk: 0.3 solar mass within 100 AU. The actual mass is somewhat lower because the disk goes isothermal after a given radius. If the temperature at 5.2 AU is 100 K, it will have fallen to 10 K at 50 AU and below the
cosmic microwave background temperature at 200 AU. As the surrounding temperatures of the molecular cloud are of
10-20 K, the 1/r fall of the temperature has to stop at some point. Boss (2002) and Pinotti et al. (2005) account for this by using a temperature profile of
,
where the parameter
TB works as a background temperature (see also
Papaloizou & Terquem 1999), leading to an essentially
isothermal profile after a certain distance. The
constancy of the mass accretion rate then dictates that the slope of
the surface density has to change accordingly, to -1.5. This, in turn, leads
to a lower disk mass. Taking into account this isothermal profile
in the very outer disk, using TB=10 K, and the corresponding change of slope in
,
the disk sets at a constant value of Q=1 after 50 AU, with a mass of 0.25 solar mass at 100 AU. If TB=20 K, the mass decreases to 0.16
at 100 AU, with Q=3 after 25 AU.
Observational studies in contrast obtain a distribution of disk masses with median of
(Andrews & Williams 2005). However, considering
an interstellar dust-to-gas ratio of 10-2, such disks contain only 15
of solid material, and cannot form the solar system. We take it as an indication that or the masses are
underestimated (Hartmann 2008) or that, as these observations correspond to relatively older disks (1-4 Myr), the gas has already dissipated significantly (Hueso & Guillot 2005). Our model then should be representative of young disks, therefore
corresponding to an early formation of relatively massive planetary embryos.
8 Summary and conclusions
We have undertaken simulations of disks of gas and solids, where the solids are represented by Lagrangian particles of radii 1, 10, 30, and 100 cm. We show that in the borders
of the dead zone, where the accretion flow stalls due a difference in turbulent viscosity, rapid gravitational collapse of the particles into protoplanets occurs within the vortices that form due to the excitation of the Rossby wave instability (Lovelace et al. 1999). As shown in LJKP08, over 300 gravitationally bound planetary embryos were formed, 20 of them being
more massive than Mars. The mass spectrum follows a power law of index
in the interval
.
Although for the main results of this study we have used sharp viscosity jumps
to model the transition between the active and dead zones, we show that the RWI
is excited up to viscosity jumps as smooth as
where H is the
pressure scale height. For this conclusion, we used the Solberg-Høiland
criterion as providing a conservative estimate of whether the RWI is excited.
The consequence of increasing the width of the viscosity jumps seems to be that
the threshold of instability is reached at increasingly longer
times. It only takes five orbits with
(the usual width used
in the models presented in this paper), but takes 350 orbits for
.
We also assessed if the vortices would survive in the more
realistic environment of a turbulent disk by making the location of the
viscosity shift oscillate with an amplitude typical of the scale length of MRI
turbulence over the period of a Keplerian revolution. We find that this has
little effect on the excitation of the RWI and saturation into vortices.
We model the solid phase with Lagrangian superparticles representing
physical pebbles and rocks of different size (1, 10, 30, and 100 cm). As these
particles are subject to different gas drag, an aerodynamical sorting by size
takes place within the vortices. The consequence of this is that the first
bound structures are formed of single particles species. This is a very
interesting result, since it is an observed fact that planetesimals are formed
of similar-sized building blocks (Scott & Krot 2005). These building blocks
seem to be sub-mm sized grains, but different nebula parameters could work as
to trap smaller particles. Youdin (2008) also points that the stationary point
of a particle trapped in vortical motion is shifted azimuthally with respect to
the eye, according to its radius .
We indeed see that clumps of
particles of different size, which collapse into different embryos inside the
same vortex, usually retain significant azimuthal displacements between each
other for long times instead of forming a single, more massive, embryo at the
vortex eye. This may or may not be a result of the size-dependent azimuthally
shifted stationary points of Youdin (2008).
A collapsed embryo is observed to be very compact. The compactness is mainly
provided by the drag force, which makes the system very dissipative (dynamical
cooling by gravity alone works on much longer timescales). Collapse towards
``infinite'' density is seen to occur in some cases, with most of the particles
occupying the same position in space (limited by numerical precision). In the
specific case when the particles dominate the time-step, the Courant condition
leads a particle to overshoot the center of the mass, so that it executes
oscillations about it, which in turn leads to a finite rms radius. We also
observe that a clump of particles is susceptible to the
disrupting effects of ram pressure erosion and gravitational tides
from the gas. Both effects are proportional to the local gas-to-solids density
ratio. When the vortices in the inner border of the dead zone undergo the transition
from the m=3 to the m=2 mode, their surface density increases, with drastic
consequences for the embryos within them. The most massive embryo by that time,
a protoplanet 6.7 times the mass of Mars, mostly formed of
cm
particles, was split into two smaller objects, of 5.9 and 0.8
,
due to the action of the gas tides. The fate of the embryos formed mostly of
10 cm was more dramatic. As the 10 cm particles experience stronger drag
forces, the ram pressure is also stronger. During the mode transition, the
combined effects of tides and erosion completely obliterated these embryos,
leaving extended arcs of particles that did not reaccumulate until the end of
the simulation. We anticipate that this effect will be very reduced in 3D
simulations. In 2D simulations, the ratio of the vertically integrated solids
density to the gas column density
never gets much
above unity even for the most massive embryo. In contrast, the ratio of the
bulk density of solids to the volume gas density
/
is greatly
increased in the midplane of 3D disks due to sedimentation.
We also observe that the solids modify the evolution of the RWI. We are drawn
to this conclusion because a simulation without the backreaction of the drag
force from the particles onto the gas developed vortices at later times when
compared to the ones that include particle feedback. We stress that this is not
due to the streaming instability, since the solids-to-gas ratio was much lower
than the value needed to excite it (
). Instead, it
is more likely that the backreaction of the drag force contributes
non-negligibly to the dispersion relation of the RWI. Self-gravity is also
seen to play a role on modifying the evolution of the turbulence. We observe
that in simulations without self-gravity, the disk tends to show less vortices
at later times. In a simulation where we switched off the self-gravity after
the vortices had developed, the dominant m=4 mode in the outer edge of the
dead zone rapidly turned into a m=3 mode. The vortices in the inner edge also
quickly turned from a dominant m=3 mode to displaying a m=2 mode instead.
We measured the Toomre Q parameter and found that the vortices have
.
This constitutes further evidence that in the presence of
self-gravity, vortex growth is not only limited by the pressure scale height
but also by the Jeans length (Mamatsashvili & Rice 2009).
An important finding in this paper is that under vortex trapping, the collisional speeds between particles are greatly reduced. We measured the collisional velocity history of every particle that is bound to the most massive embryo at the end of the simulation, and found that the vast majority of them never experienced close encounters at speeds greater than 1 m s-1. This is well below the fragmentation threshold and lends support to the long-held idea that vortices provide a superbly favorable environment for planetary growth. Growth by coagulation beyond the optimal size for planet formation is also avoided because the timescales for coagulation are much longer than the rapid timescale for gravitational collapse witnessed in our models. This does not exclude the possibility that coagulation itself is sped up within a vortex. In this case, the vortex will behave as a planetesimal factory, quickly producing kilometer sized bodies that leave the vortex due to their decoupling from the gas. This, as noted by Klahr & Bodenheimer (2006) is very positive for planet formation, since the formed planetesimals are then scattered through the disk, where they can be used to form planets independently of a vortex. Even though it implies that we are facing the comfortable position of a win-win situation for planet formation, one has to decide which process (planet formation by direct gravitational collapse or planetesimal formation by fast coagulation) is getting the upper hand inside the Rossby vortices. A definite answer to this question can only be drawn from a simulation that includes the processes of coagulation/fragmentation. Unfortunately, inclusion of sophisticated coagulation/fragmentation models such as that of Brauer et al. (2007) would render a hydrodynamical simulation overly expensive. A possible alternative would be a Monte Carlo model of dust coagulation, such as the one recently developed by Ormel et al. (2007). A further development of the Monte Carlo technique is described in Zsom & Dullemond (2008). The main difference between the two models is that while Ormel et al. (2007) simulate coagulation between real dust particles, Zsom & Dullemond (2008) use superparticles to model coagulation and fragmentation. Therefore the latter one is more suitable for hydrodynamical simulations such as ours and simple estimations show that this model could be adapted to a hydro model with no prohibitive overhead.
We reiterate that the models presented suffer from a number of limitations, detailed in Sect. 7. Some of them, like refining the particle mass resolution to the individual pebbles and rocks, are beyond the capabilities of the current generation of computer models. Others, however, such as inclusion of detailed thermal physics, could be tackled with relatively little effort. We urge researchers active on the field to consider these problems. It is our hope that a coherent picture of planet formation in the magnetically dead zones of accretion disks shall emerge as a result of it.
Acknowledgements
Simulations were performed at the PIA cluster of the Max-Planck-Institut für Astronomie and on the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX). This research has been supported in part by the Deutsche Forschungsgemeinschaft DFG through grant DFG Forschergruppe 759 ``The Formation of Planets: The Critical First Growth Phase''. A. Zsom acknowledges support by the IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg.
Appendix A: Resolution study
In Fig. A.1 we present a resolution study of our models. The upper panels show Cartesian models, including self-gravity and solids, while the lower panels show cylindrical models with dustless non-selfgravitating gas. Both are shown in selected snapshots at 10, 20, and 30 orbits.
In the Cartesian runs, the vortices in the high-resolution run (512)
are excited earlier than in the middle resolution run (256
).
At later times, the vortices in the high-resolution run also
appear sharper. While in the two runs the vortices in the outer edge
of the dead zone look remarkably similar, the inner edge displays
very different behavior. Up to 30 orbits the middle resolution run has not
shown signs of prominent vortices. In contrast, the high resolution run
had the
inner edge developing vortices as early as 10 orbits. This is not only
due to the high resolution run having a wider inertial range, but also
because a flow with cylindrical symmetry is more coarsely resolved
near the origin in a Cartesian grid. Because of this, the most unstable
modes are under-resolved in the inner disk of the middle resolution run.
A Cartesian run with low resolution (128
,
not shown) did not
develop vortices even in the outer disk by the same time. At 30 orbits,
the density inside the vortices peak at similar values,
and 3.7 for the middle and high resolution runs, respectively.
The cylindrical runs also illustrate the small amount of differences between
the vortices in different runs. With better azimuthal resolution, the
vortices are excited even in the 128
run, and in both runs they
peak with surface density
.
In fact, the main effect
of resolution appears to be that, as it increases, the RWI is excited
increasingly earlier. The high-resolution run displays weaker vortices
than the others because in this case we were forced to use shock viscosity.
![]() |
Figure A.1:
The development of the dead zone vortices in different resolution and grid geometry. With increased resolution the RWI is excited increasingly earlier. In general the vortices in cylindrical runs look sharper than in the Cartesian ones, due to the better azimuthal resolution. The Cartesian run of middle resolution (256 |
Open with DEXTER |
Appendix B: Emulating turbulent motions
![]() |
Figure B.1: Evolution of the RWI with a time-varying viscosity jump as specified by Eq. (13) and Eqs. (B.1)-(B.2). The panel on the left-hand side measures the violation of the Solberg-Høiland criterion. The right-hand side panel shows the appearance of the disk at 30 orbits. Both the inner and outer edge quickly reach the threshold of instability ( left panel). At 30 orbits, the inner edge already reached a saturated state and launched vortices ( right panel). |
Open with DEXTER |
In this study, we considered the dead zone to be represented by a static viscosity profile. In a more realistic scenario, turbulent motions caused by the MRI and variations in the coupling between the magnetic field and the plasma will give rise to a turbulent resistivity. This is expected to cause the border of the dead zone to vary in space and time, with implications for the evolution of the RWI.
To assess the impact of space and time variability of the edges of the dead zone, we model the viscous jumps using
but make the jumps shift radially in time by substituting r1 and r2 in Eq. (13) by
where h=H/r is the aspect ratio. So, the location shifts by two scale heights over a Keplerian revolution. The results are shown in Fig. B.1, where we show the appearance of the disk at 30 orbits and the azimuthal average of

In this simulation, the RWI is still excited and vortices are launched. The main difference when compared to simulations with static profiles is that the instability takes more time
to violate the Solberg-Høiland criterion, 10 orbits, compared to 5 in the static case. This is due to the fact that the shifting viscosity jump smears the pressure maximum, so the amplitude of the pressure jump is shorter and the width is larger than in the static case.
Appendix C: Onset of erosion
![]() |
Figure C.1:
A clump of 1000 massive particles moving against headwinds
of
|
Open with DEXTER |
According to Eq. (30), a tight distribution of particles under Epstein drag
should suffer erosion if
![]() |
(C.1) |
In this appendix, we perform numerical simulations to test the validity of this condition. We model a clump of 103 particles suffering a strong headwind from the gas (




We plot in Fig. C.1 the evolution of clumps for four different values of
.
For
and
,
is less than 1 so the
clump is stable against ram pressure and contracts. The other clumps (
and
)
have
above unit, and experience intense erosion.
We also considered a flow of Mach number
.
In this case the initial
gravitational Weber number is
,
and the
clumps of smaller particles are supposed to be more stable. Indeed, this is what we
see in the figure. The clump of particles of
is now marginally
stable and contracts, experiencing much less erosion than in the
case.
We would like to draw attention to an interesting feature of Fig. C.1.
The cases of
and
(upper panels)
provide yet another perspective for the action of drag force cooling. Both clumps are
stable against erosion but the clump of
has shrunk
considerably more than the clump of
,
which looks very
extended. What is happening is that the clump with
is
too weakly coupled to the gas and therefore takes longer to collapse, as
seen in the time series of the rms position (Fig. C.1, lower panels).
References
- Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134 [NASA ADS] [CrossRef] (In the text)
- Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 [NASA ADS] [CrossRef] (In the text)
- Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1 [NASA ADS] [CrossRef] (In the text)
- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] (In the text)
- Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [NASA ADS] [CrossRef] (In the text)
- Barge, P., & Sommeria, J. 1995, A&A, 295, 1 [NASA ADS] (In the text)
- Benz, W. 2000, SSRv, 92, 279 [NASA ADS] (In the text)
- Boss, A. P. 2002, ApJ, 576, 462 [NASA ADS] [CrossRef] (In the text)
- Bracco, A., Chavanis, P.-H., Provenzale, A., & Spiegel, E. A. 1999, Phys. Fluids, 11, 2280 [NASA ADS] [CrossRef] (In the text)
- Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Brauer, F., Dullemond, C. P., & Henning, Th. 2008a, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Brauer, F., Henning, Th., & Dullemond, C. P. 2008b, A&A, 487, 1 [NASA ADS] [CrossRef] [EDP Sciences]
- Carballido, A., Stone, J. M., & Turner, N. J. 2008, MNRAS, 386, 145 [NASA ADS] [CrossRef] (In the text)
- Chavanis, P.-H. 2000, A&A, 356, 1089 [NASA ADS] (In the text)
- Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [NASA ADS] [CrossRef] (In the text)
- Crow, S. C. 1970, AIAAJ, 8, 2172 [NASA ADS] [CrossRef] (In the text)
- Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432 [NASA ADS] [CrossRef] (In the text)
- Dominik, C., Blum, J., Cuzzi, J. N., & Wurm, G. 2007, Protostars and Planets V, 783 (In the text)
- Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] (In the text)
- Edgar, R. G., & Quillen, A. C. 2008, MNRAS, 387, 387 [NASA ADS] [CrossRef] (In the text)
- Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 [NASA ADS] [CrossRef] (In the text)
- Fromang, S., & Nelson, R. P. 2005, MNRAS, 364, 81 [NASA ADS] (In the text)
- Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
- Garaud, P., & Lin, D. N. C. 2004, ApJ, 608, 1050 [NASA ADS] [CrossRef] (In the text)
- Gledzer, E. B., Dolzhanskii, F. V., Obukhov, A. M., & Ponomarev, V. M. 1975, Isv. Atmos. Ocean. Phys., 11, 617
- Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] (In the text)
- Haghighipour, N., & Boss, A. P. 2003, ApJ, 583, 996 [NASA ADS] [CrossRef] (In the text)
- Hartmann, L. 2008, Phys. Scrip., T130, 4012 [NASA ADS] [CrossRef] (In the text)
- Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] (In the text)
- Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Johansen, A., Henning, Th., & Klahr, H. 2006, ApJ, 643, 1219 [NASA ADS] [CrossRef] (In the text)
- Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] (In the text)
- Johansen, A., Brauer, F., Dullemond, C., Klahr, H., & Henning, T. 2008, A&A, 486, 597 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Johnson, B. M., & Gammie, C. F. 2005, ApJ, 635, 149 [NASA ADS] [CrossRef] (In the text)
- Kerswell, R. R. 2002, Annu. Rev. Fluid Mech., 34, 83 [NASA ADS] [CrossRef] (In the text)
- Klahr, H. 2004, ApJ, 606, 1070 [NASA ADS] [CrossRef] (In the text)
- Klahr, H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] (In the text)
- Klahr, H., & Bodenheimer, P. 2006, ApJ, 639, 432 [NASA ADS] [CrossRef] (In the text)
- Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] (In the text)
- Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874 [NASA ADS] [CrossRef] (In the text)
- Lodato, G. 2008, New Astron. Rev., 52, 21 [NASA ADS] [CrossRef] (In the text)
- Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] (In the text)
- Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008a, A&A, 479, 883 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008b, A&A, 491, L41 [NASA ADS] [CrossRef] [EDP Sciences] (LJKP08) (In the text)
- Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Lyttleton, R. A. 1972, MNRAS, 158, 463 [NASA ADS] (In the text)
- Mamatsashvili, G. R., & Rice, W. K. M. 2009 [arXiv:0901.1617] (In the text)
- Marcus, P. 1990, Journal of Fluid Mechanics, 215, 393 [NASA ADS] [CrossRef] (In the text)
- Méheut, H., Casse, F., Varnière, P., & Tagger, M. 2008, Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Charbonnel, F. Combes, & R. Samadi, 417 (In the text)
- Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398 [NASA ADS] [CrossRef] (In the text)
- Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805 [NASA ADS] [CrossRef] (In the text)
- Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences]
- Paardekooper, S.-J. 2006, PhDT (In the text)
- Papaloizou, J., & Terquem, C. 1999, ApJ, 521, 823 [NASA ADS] [CrossRef] (In the text)
- Paraskov, G. B., Wurm, G., & Krauss, O. 2006, ApJ, 648, 1219 [NASA ADS] [CrossRef] (In the text)
- Petersen, M. R., Julien, K., & Stewart, G. R. 2007, ApJ, 658, 1236 [NASA ADS] [CrossRef] (In the text)
- Pinotti, R., Arany-Prado, L., Lyra, W., & Porto de Mello, G. F. 2005, MNRAS, 364, 29 [NASA ADS] [CrossRef] (In the text)
- Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543 [NASA ADS] [CrossRef] (In the text)
- Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2006, MNRAS, 372, 9 [NASA ADS] (In the text)
- Safronov, V. S. 1969, QB, 981, 26 (In the text)
- Scott, E. R. D., & Krot, A. N. 2005, in ASP Conf. Ser., 341, 15 (In the text)
- Sekiya, M. 1998, Icarus, 133, 298 [NASA ADS] [CrossRef] (In the text)
- Shakura, N. I., & Suyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
- Shen, Y., Stone, J. M., & Gardiner, T. A. 2006, ApJ, 653, 513 [NASA ADS] [CrossRef] (In the text)
- Sirono, S.-I. 2004, Icarus, 167, 431 [NASA ADS] [CrossRef] (In the text)
- Tanga, P., Babiano, A., Dubrulle, B., & Provenzale, A. 1996, Icarus, 121, 158 [NASA ADS] [CrossRef] (In the text)
- Terquem, C. 2008, ApJ, 689, 532 [NASA ADS] [CrossRef] (In the text)
- Turner, N. J., & Sano, T. 2008, ApJ, 679, 131 [NASA ADS] [CrossRef]
- Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69 [NASA ADS] [CrossRef]
- de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] (In the text)
- de Val-Borro, M., Artymowicz, P., D'Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043 [NASA ADS] [CrossRef] [EDP Sciences]
- Varnière, P., & Tagger, M. 2006, A&A, 446, 13 [NASA ADS] [CrossRef] (In the text)
- Weidenschilling, S. J. 1977a, MNRAS, 180, 57 [NASA ADS] (In the text)
- Weidenschilling, S. J. 1977b, Ap&SS, 51, 153 [NASA ADS] [CrossRef]
- Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] (In the text)
- Weidenschilling, S. J., & Cuzzi, J. N. 1993, Protostars and Planets III, 1031 (In the text)
- Woitke, P., & Helling, C. 2003, A&A, 399, 297 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
- Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] (In the text)
- Youdin, A. N. 2008 [arXiv:0807.1114] (In the text)
- Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] (In the text)
- Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] (In the text)
- Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
Footnotes
- ...
- See http://www.nordita.org/software/pencil-code
- ... transport
- Another
alternative is local ionization provided by the decay of the short-lived
radioactive nuclide
Al (Umebayashi & Nakano 2009), although Turner & Sano (2008) show that the free electrons given out by this low ionization source would quickly recombine on the surface of
m-sized dust grains.
- ...
experience
- The Reynolds number of the flow past a clump can be written as
. At the grid scale our choice of viscosity is usually
cm s
(it decreases very fast as we go to smaller wavenumbers, as k6). For a clump as the ones of this study, of
km and
m s
, the Reynolds number is
. At such incredibly low Reynolds numbers, inertia plays no role. The self-gravity of the particles, therefore, is not holding the cluster together against drag forces from the gas, but against largely exaggerated viscous stresses. In comparison, in the solar nebula, the molecular viscosity is much lower and the Reynolds number is expected to be >106 (Cuzzi et al. 2008).
All Figures
![]() |
Figure 1:
Enlargement around one of the vortices in a snapshot at 20
orbits. Contours of |
Open with DEXTER | |
In the text |
![]() |
Figure 2: Evolution of the pressure bump with different widths of the viscosity jump ( upper most panels; the dashed line is the Heaviside jump, for comparison). The violation of the Solberg-Høiland criterion ( lower most panels) is a conservative indication that the threshold of instability of the RWI was reached. As the width of the jump increases, the threshold takes increasingly longer time to be breached. For jumps smoother than two scale lengths, the threshold was not reached (up to 500 orbits). |
Open with DEXTER | |
In the text |
![]() |
Figure 3: 1D simulation of collapse of 1000 particles. With drag force the kinetic energy of the particles is efficiently dissipated and the particles collapse at subgrid scale towards infinite density. When the drag force is excluded the system cannot dissipate energy and a halo of particles, 10 grid cells wide, is formed. |
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Evolution of the most massive clump formed (solid line), which we traced back in time from the end of the simulation. It differs from the instantaneous most massive clump (dashed line) because the clumps have different feeding rate and can also experience mass loss, as in the episode that happened at |
Open with DEXTER | |
In the text |
![]() |
Figure 5: Time series of the mass loss episode. Due to gravitational tides from the massive vortex whence the embryo formed, a large chunk of particles was detached from the original cluster. At 87 orbits, a separation 4 times the Earth-Moon distance is seen. The separation grows and two orbits later the two bodies do not overlap Hill radii, thus counting as separate objects. |
Open with DEXTER | |
In the text |
![]() |
Figure 6: The state of the disk before a)- d) and after e)- h) the mass loss episode. The conspicuous difference between them is that the embryo has left its parental vortex from one snapshot to the other. It is seen as a bright spot in panels f) and h), at (X,Y)=(-0.65, -0.19). In panel b) (before the mass loss), the embryo is at (X,Y)=(-0.40, -0.53) but not easily spotted among the swarm of solids inside the vortex. Panel i) shows a horizontal slice through this location, in which we see that the density of solids does not peak much higher than the gas density at the location of the embryo (panel j)). Significant gas tides are expected as the gravitational potential (panel k)) and acceleration (panel l)) have similar contributions from the gas and solid components. |
Open with DEXTER | |
In the text |
![]() |
Figure 7:
Can self-gravity alone maintain the formed embryos together? We switch off the hydrodynamics at the last snapshot and run for additional 200 orbits. The most massive embryo, which had already left the parental vortex (at |
Open with DEXTER | |
In the text |
![]() |
Figure 8:
The relative velocities between particles are contrasted with
collisional velocities calculated from Eq. (32). At the
initial condition (t=0) the relative velocities between particles follow
what is expected from the Keplerian shear (dashed line), reaching as far
as 300 m s-1 inside a grid cell (which is 31
|
Open with DEXTER | |
In the text |
![]() |
Figure 9:
a) Accretion history of the most massive embryo.
Collapse happens at 40 orbits, and further accretion of particles happens through the next 18 orbits, after which the maximum mass is attained. This plot defines the time t0 that each particle is accreted. b) Histograms of the maximum collisional velocities experienced by a given particle, before (dot-dashed line) and after (solid line) its accretion time t0. The solid line histogram represents the maximum collisional speed from t0 until the end of the simulation. The dot-dashed line histogram covers a time interval of ten orbits before t0. As the vortices form at |
Open with DEXTER | |
In the text |
![]() |
Figure 10:
Collisional speeds as a function of separation at selected snapshots at 30, 35, 40, and 45 orbits (
|
Open with DEXTER | |
In the text |
![]() |
Figure 11: Aerodynamical sorting for the particles trapped in the vortical motion. The figure is centered at the vortex shown in Fig. 1, at 18 orbits. The unit of length is the Earth-Moon mean distance, and the Y-coordinate points to the star. As particles of different size have different friction times, differential drag occurs inside the vortex, effectively sorting the particles spatially by size. |
Open with DEXTER | |
In the text |
![]() |
Figure 12:
The radius of the particle subject to maximum drift (
|
Open with DEXTER | |
In the text |
![]() |
Figure 13: Time series of the gas and solid density (for each individual particle size), showing the destruction of the embryos of 10 cm particles when the inner vortices undergo the transition from the m=3 to the m=2 mode. At 75 orbits ( upper panels), embryos are seen in both the inner and outer vortices, for both 10 cm and 30 cm particles. At 85 orbits one of the embryos in the 10 cm phase was disrupted, followed by a second at 95 orbits, and the last one ten orbits later. The embryos composed of 30 cm particles also experience tides. But as their gravitational Weber numbers are smaller, they just undergo splitting, the large fragments being more stable against erosion than the embryos composed of particles of 10 cm. |
Open with DEXTER | |
In the text |
![]() |
Figure 14: Upper panels: evolution of a disk without particles and without self-gravity, which serves as a control run for the next plots. Middle panels: evolution of a disk without particles but with self-gravity. The difference compared to the upper panels is that the dominant mode in the outer disk changed from m=2 to m=5. Self-gravity modifies the dispersion relation of the RWI, or it stalls the inverse cascade of power known to occur in 2D turbulence, or both. Lower panels: evolution of a disk without self-gravity but with particles. The backreaction leads to an early excitation of the RWI in the inner edge of the dead zone. Conversely, the outer edge goes unstable later when compared to the other two runs. Since the particle density is not high enough to excite the streaming instability, we take it as evidence that the backreaction modifies the dispersion relation of the RWI. |
Open with DEXTER | |
In the text |
![]() |
Figure 15: In the simulation shown in LJKP08, the vortices maintained a m=4 mode in the outer edge of the dead zone until the end of the simulation. By switching self-gravity off, we see that in less than 15 orbits the outer vortices merged into a m=3 mode, while the m=3 mode in the inner disk edge into a m=2 mode. This is because as Q decreases, the size of the vortices approaches the Jeans length scale, which effectively halts the inverse cascade of energy. |
Open with DEXTER | |
In the text |
![]() |
Figure A.1:
The development of the dead zone vortices in different resolution and grid geometry. With increased resolution the RWI is excited increasingly earlier. In general the vortices in cylindrical runs look sharper than in the Cartesian ones, due to the better azimuthal resolution. The Cartesian run of middle resolution (256 |
Open with DEXTER | |
In the text |
![]() |
Figure B.1: Evolution of the RWI with a time-varying viscosity jump as specified by Eq. (13) and Eqs. (B.1)-(B.2). The panel on the left-hand side measures the violation of the Solberg-Høiland criterion. The right-hand side panel shows the appearance of the disk at 30 orbits. Both the inner and outer edge quickly reach the threshold of instability ( left panel). At 30 orbits, the inner edge already reached a saturated state and launched vortices ( right panel). |
Open with DEXTER | |
In the text |
![]() |
Figure C.1:
A clump of 1000 massive particles moving against headwinds
of
|
Open with DEXTER | |
In the text |
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.