---

Planet formation bursts at the borders of the dead zone in 2D numerical simulations of circumstellar disks

W. Lyra¹ - A. Johansen² - A. Zsom³ - H. Klahr³ - N. Piskunov¹

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, and are the vertically integrated gas density and bulk density of solids, respectively; stands for the velocity of the gas parcels; is the position and is the velocity of the solid particles, P is the vertically integrated pressure, is the sound speed, the gravitational potential, the viscosity, and the rate-of-strain tensor. The functions and are sixth order hyperdiffusion and hyperviscosity terms that provide extra dissipation near the grid scale, explained in Lyra et al. (2008a). They are needed because the high order scheme of the Pencil Code has too little overall numerical dissipation.

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. C_D 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 and are the coefficients of Epstein and Stokes drag, respectively. These are

where is the Mach number, is the Reynolds number of the flow past the particle, and is the kinematic viscosity of the gas. We defer the reader to Lyra et al. (2009) for further details of the drag force. The self-gravity solver is also explained in that paper.

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 . We choose c_s0=0.05 and r₀=5.2 AU. We impose that the disk should have 30   of solid material within the modeled range. We also use a surface density of solids following the same slope of the surface density , where is the initial solids-to-gas ratio. For the interstellar ratio , the surface density at r₀ is  g cm (which we round to  g cm). The gas disk thus has a mass of within the modeled range. We further discuss the disk mass in Sect. 7.5. The gas has a Toomre Q parameter of 30 at r₀, and is throughout stable against gravitational instability.

The dead zone is modeled as static viscosity jumps following arc-tangent profiles

where r₁=0.6 and r₂=1.2 are the locations of the jumps and its width. We adopt , which provides a smooth jump over two grid cells only. The jump is thus close to a Heaviside function yet still differentiable. We use in code units, corresponding to an alpha value (Shakura & Sunyaev 1973) of .

For the solids, we use 10⁵ 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

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

The quantity is defined as

where

is the epicyclic frequency and

are the radial length scale of the entropy and density variations, respectively. is the adiabatic index. For corotational modes ( ) in a barotropic ( ) disk, the quantity reduces to , which is readily identified with (half) the inverse of vortensity 

which in turn led Lovelace et al. (1999) to interpret as an entropy-modified version of, or generalized, potential vorticity. An extremum in the profile of can be generated, for example, by a pressure bump somewhere in the disk. The dispersion relation of the disturbances is analogous to the dispersion relation of Rossby waves in planetary atmospheres, hence the name of the instability.

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 suffices to determine stability. The other term

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 N² 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), is the vertically integrated specific entropy of the gas. The function is a sixth order hyperconductivity term, analogous to hyperdiffusion for density. The adiabatic index is .

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 r₀. 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

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, is the vector radius of each particle and is the vector radius of their center of mass. We take this value as a measurement of the ``radius'' of the embryo. The most massive embryo has a radius of 1.13 . This compactness corresponds to a tenth of a thousandth of the resolution element.

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.

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.

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.

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

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 F_G = -GM/R² for the planet's (self-)gravitational force at its surface. Substituting , we have , so

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 / is expected to be much lower than the ratio of the column gas density to the vertically integrated surface 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, and the exaggerated viscous stresses might de-stabilize potentially stable clumps.

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

For Epstein drag (Eq. (10)), C_D does not depend on . So, for all other quantities being constant, we expect to decrease linearly as increases. In other words, self-gravitating clumps of larger particles should be more stable than clumps composed of smaller particles.

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 , a clump of at will be stable if . The embryos we consider are formed predominantly of particles of 10 cm ( ) and 30 cm ( ). We conclude that erosion might play a role in our case. The sharp dependence of on the distance in Eq. (31) also means that embryos at the inner edge of the dead zone are more prone to erosion than the ones at the outer edge. We will develop this further in Sect. 5.2.2.

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 r_ij between them tends to zero. For our resolution and choice of r₀, 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 r_ij 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.

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 is much smaller than the shear dependence of (where is the pressure-correction factor), we can assume that it also corrects for this quantity. Figure 8 also shows the collisional speeds calculated by Eq. (32). The dependence on separation was greatly suppressed.

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 t₀ 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 t₀ and the end of the simulation (200 orbits). The former refers to a time interval of 10 orbits before t₀. 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.

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

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.

5.2 Differences between the inner and outer embryos

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

and calculating the radius for a given stopping time . Substituting (Eq. (10) at subsonic regime) and Stokes law (with ) yields the quadratic equation for 

and, as the radius is positive, the solution is unique

Figure 12 shows this curve as a function of radius for and our initial parameters. The particle radius of predicted by both pure Epstein and pure Stokes drag are shown for comparison. The mean free path  is also shown. The figure shows that the curve is quite flat for our choice of parameters, so inner and outer vortices must have similar efficiency on trapping particles of a given size.

5.2.2 Tidal disruption and erosion of the inner embryos of 10 cm

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).

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 T_B 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 T_B=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 T_B=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.

Appendix B: Emulating turbulent motions

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 r₁ and r₂ 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 the same 2D model, as opposed to 1D as in Sect. 3.2).

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

According to Eq. (30), a tight distribution of particles under Epstein drag should suffer erosion if

In this appendix, we perform numerical simulations to test the validity of this condition. We model a clump of 10³ particles suffering a strong headwind from the gas ( ). The blob of particles is initially set as a tight Gaussian distribution about the center of the grid, with surface density peak of . In units where , the initial gravitational Weber number at the surface of the blob is therefore .

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 k⁶). 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 >10⁶ (Cuzzi et al. 2008).

All Figures

---

Copyright ESO 2009