On the dynamics of planetesimals embedded in turbulent protoplanetary discs with dead zones
aa r X i v : . [ a s t r o - ph . E P ] J un Mon. Not. R. Astron. Soc. , 1–18 (2002) Printed 22 June 2018 (MN L A TEX style file v2.2)
On the dynamics of planetesimals embedded in turbulentprotoplanetary discs with dead zones
Oliver Gressel ⋆ , Richard P. Nelson ⋆ and Neal J. Turner ⋆ Astronomy Unit, Queen Mary, University of London, Mile End Road, London E1 4NS, United Kingdom Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
Accepted 2011 April 19. Received 2011 April 19; in original form 2011 February 28
ABSTRACT
Accretion in protoplanetary discs is thought to be driven by magnetohydrodynamic (MHD)turbulence via the magnetorotational instability (MRI). Recent work has shown that a plan-etesimal swarm embedded in a fully turbulent disc is subject to strong excitation of the veloc-ity dispersion, leading to collisional destruction of bodies with radii R p <
100 km. Significantdi ff usion of planetesimal semimajor axes also arises, leading to large-scale spreading of theplanetesimal population throughout the inner regions of the protoplanetary disc, in apparentcontradiction of constraints provided by the distribution of asteroids within the asteroid belt.In this paper, we examine the dynamics of planetesimals embedded in vertically stratifiedturbulent discs, with and without dead zones. Our main aims are to examine the turbulentexcitation of the velocity dispersion, and the radial di ff usion, of planetesimals in these discs.We employ three dimensional MHD simulations using the shearing box approximation, alongwith an equilibrium chemistry model that is used to calculate the ionisation fraction of the discgas as a function of time and position. Ionisation is assumed to arise because of stellar X-rays,galactic cosmic rays and radioactive nuclei. In agreement with our previous study, we find thatplanetesimals in fully turbulent discs develop large random velocities that will lead to colli-sional destruction / erosion for bodies with sizes below 100 km, and undergo radial di ff usionon a scale ∼ . R p
100 km. Wealso find that radial di ff usion occurs over a much reduced length scale ∼ .
25 au over the disclife time, this being consistent with solar system constraints. We conclude that planetesimalgrowth via mutual collisions between smaller bodies cannot occur in a fully turbulent disc.By contrast, a dead zone may provide a safe haven in which km-sized planetesimals can avoidmutual destruction through collisions.
Key words: accretion disks – magnetohydrodynamics (MHD) – methods: numerical – plan-etary systems: formation – planetary systems: protoplanetary disks
The formation of planetesimals is a key stage in the formation ofplanetary systems, but at present there is little consensus aboutthe processes involved. One class of models envisages an incre-mental process in which small dust grains collide, stick and set-tle vertically within a protoplanetary disc, with continued growththrough coagulation eventually forming km-sized planetesimals(Weidenschilling & Cuzzi 1993; Weidenschilling 2000). Modelsdiscussed by Weidenschilling (2000), that assume the disc is lami-nar and has a surface density a few times larger than the minimummass solar nebula model of Hayashi (1981), predict planetesimalformation times that are a few thousand local orbital periods at each ⋆ E-mail: { o.gressel, r.p.nelson } @qmul.ac.uk; [email protected] location in the disc, implying formation times of ≈ few × yr at5 au. More recent models presented by Brauer et al. (2008), that ex-amine planetesimal formation in low density turbulent discs, con-clude that incremental planetesimal formation is prevented by therapid inward migration of metre-sized boulders when they form,combined with fragmentation induced by high velocity collisionscaused by di ff erential radial migration. Large km-sized bodies wereonly obtained in these models when unrealistically large values forthe critical velocity for fragmentation were adopted, and a signif-icant increase in the dust-to-gas ratio above solar abundance wasassumed. A local particle-in-a-box model of growth through dustcoagulation presented by Zsom et al. (2010), that incorporates re-cent experimental results on the compactification of porous aggre-gates, suggests that coagulative growth beyond millimetre sizes isdi ffi cult to achieve because compacted aggregates tend to bounce c (cid:13) Gressel, Nelson & Turner rather than stick. These latter results raise serious questions aboutthe validity of the incremental picture of planetesimal formation.A number of the di ffi culties faced by the incremental modelhave been known for many years, leading to searches for alterna-tive planetesimal formation scenarios. Gravitational instability ofthe dense dust layer that forms through vertical settling was sug-gested by Goldreich & Ward (1973) as a means of forming km-sized planetesimals. Although this pathway is not now widely ac-cepted, recent models of planetesimal formation that involve col-lective e ff ects, including self-gravity of regions of enhanced solidsconcentration, have been proposed. Cuzzi et al. (2008) have de-veloped a model in which mm-sized chondrules are concentratedin turbulent eddies, and then contract on the settling time of thesolid grains under the action of self-gravity to form large 100 km-sized planetesimals. The combined action of the streaming instabil-ity (Youdin & Goodman 2005), and trapping of metre-sized boul-ders in short-lived turbulent eddies, has been suggested as a meansof inducing planetesimal formation via gravitational collapse byJohansen et al. (2007), leading to the direct formation of planetesi-mals that are as large as the asteroid Ceres ( R p ≃
500 km). Althoughthese methods of forming planetesimals avoid the problem of grow-ing beyond the metre-sized barrier, they each have their own di ffi -culties, and as yet no clear consensus of how planetesimals formhas emerged.Once planetesimals have formed, runaway growth may leadto the formation of ∼ ∼ yr(Wetherill & Stewart 1993; Kenyon & Bromley 2009), which thenaccumulate smaller bodies to become planetary embryos and coresduring the oligarchic phase (Ida & Makino 1993; Kokubo & Ida1998). Rapid runaway growth of km-sized planetesimals requiresthat their velocity dispersion remains significantly smaller than theescape velocity of the accreting bodies, which for 10 km-sized ob-jects is ∼
10 m s − . A self-stirring swarm of planetesimals embed-ded in a laminar disc can maintain a low velocity dispersion viagas drag and dynamical friction, but the situation in a turbulentdisc may be very di ff erent due to stochastic forcing of planetesimalrandom motions by the turbulence. Planetesimal accretion that oc-curs at a rate determined by the geometric cross section may leadto planetary formation time scales in excess of observed disc lifetimes of a few Myr (Haisch et al. 2001).A further constraint is that planetesimal collision velocitiesmust remain below the catastrophic disruption threshold for col-lisional growth to occur. For planetesimals between the sizes 1- 10 km, disruption thresholds lie approximately in the interval10 m s − v crit
60 m s − , depending on the material strengthand whether the planetesimals are monolithic bodies or rubble piles(Benz & Asphaug 1999; Stewart & Leinhardt 2009). Clearly an ex-ternal source of stirring, such as disc turbulence, may not only pre-vent runaway growth, but may lead to the collisional destruction ofplanetesimals.Turbulence is generally believed to provide the anomaloussource of viscosity required to explain the typical ∼ − M ⊙ yr − accretion rates of T Tauri stars (Sicilia-Aguilar et al. 2004). Themost promising source of this is MHD turbulence driven by themagnetorotational instability (MRI, Balbus & Hawley 1991). Bothlocal shearing box and global simulations indicate that the non lin-ear, saturated state of the MRI is vigorous turbulence that is ableto transport angular momentum at a rate that can explain obser-vations of T Tauri accretion rates (Hawley et al. 1995; Armitage1998; Hawley 2001; Papaloizou & Nelson 2003). But the MRI re-quires that the gas is su ffi ciently ionised in order for the linear in-stability to operate (Blaes & Balbus 1994), and for non linear tur- bulence to be sustained (Fleming et al. 2000). In the absence of suf-ficient free electrons in the gas phase, the flow returns to a laminarstate.Protoplanetary discs surrounding T Tauri stars are cold anddense, leading to low levels of ionisation (Umebayashi 1983). Inthe main body of these discs the primary sources of ionisation areexpected to be external: stellar X-rays (Glassgold et al. 1997) andpossibly galactic cosmic rays (Umebayashi & Nakano 1981), eachof which have a limited penetration depth into the disc. Gammie(1996) presented a model in which the disc surface layers areionised, and sustain a turbulent accretion flow, whereas the disc in-terior is a ‘dead zone’ where the flow remains laminar and minimalaccretion takes place. This layered accretion disc model has beenthe subject of numerous investigations over recent years which haveexamined the role of dust grains (Sano et al. 2000; Ilgner & Nelson2006a), gas-phase heavy metals (Fromang et al. 2002), di ff erentchemical reaction networks (Ilgner & Nelson 2006a), the role ofturbulent mixing (Ilgner & Nelson 2006b), and the Hall e ff ect(Wardle 1999).MHD simulations have also been used to examine the dynam-ics of magnetised discs in which a finite electrical conductivityplays an important role. Fleming et al. (2000) examined the satu-rated level of turbulence as a function of the magnetic Reynoldsnumber. Fleming & Stone (2003) examined disc models in whichmagnetic resistivity varied with height, simulated the layered accre-tion model proposed by Gammie (1996), and demonstrated that thedead zone maintains a modest Reynolds stress due to waves beingexcited by the active layers. More recent studies have coupled timedependent chemical networks to the MHD evolution using a multi-fluid approach, and have demonstrated that turbulent di ff usion ofcharge carriers into the dead zone can enliven it, in the absence ofsmall dust grains (Turner et al. 2007; Ilgner & Nelson 2008). Thepresence of a significant population of small grains, however, re-duces the time for removal of free charges from the gas phase, andthe dead zone is maintained (Turner & Sano 2008).Low mass planets and planetesimals embedded in turbulentprotoplanetary discs are subject to stochastic gravitational forcesdue to turbulent density fluctuations (Nelson & Papaloizou 2004;Nelson 2005). The stochastic forcing leads to di ff usion of the semi-major axes and excitation of the eccentricity of the embedded bod-ies. Nelson & Gressel (2010) – hereafter paper I – examined thedynamics of planetesimals embedded in fully turbulent cylindricaldisc models without dead zones, and demonstrated that the exci-tation of the velocity dispersion will lead to collisional disruptionof planetesimals of size R p
10 km. They also demonstrated thatvery good agreement can be obtained between local shearing boxsimulations and global simulations when shearing boxes of su ffi -cient size are employed (see also Yang et al. 2009). In this paper,we extend the work presented in paper I and examine the dynam-ics of planetesimals embedded in vertically stratified discs, withand without dead zones, using shearing box simulations. A key re-sult is that the turbulent stirring of planetesimals is significantlyreduced in discs with dead zones, possibly allowing for their con-tinued growth rather than destruction or erosion during collisions.As such, we propose that dead zones may provide safe havens forplanetesimals, and are a required ingredient for the formation ofplanetary systems.This paper is organised as follows: We describe the physicalmodel and numerical methods in Section 2. The presentation ofthe results is divided into four parts: In Section 3, we describe thegeneral morphology and dynamics of the emerging hydromagneticturbulence. Section 4 discusses how the resulting stochastic grav- c (cid:13) , 1–18 ynamics of planetesimals in dead zones itational torques experienced by bodies embedded in the disc canbe classified by their distribution and auto-correlation. Finally, inSections 5 and 6, we analyse the temporal evolution in eccentricityand radial di ff usion of a swarm of particles immersed in the flow.We discuss the implications of our results for planetary formationin Section 7, and we draw our conclusions in Section 8. For the models reported-on in this paper, we make use of thesecond-order Godunov code nirvana - iii (Ziegler 2004, 2008) andsolve the standard equations of resistive magneto-hydrodynamics(MHD) in a local, corotating Cartesian frame (ˆ x , ˆ y , ˆ z ), which read: ∂ t ρ + ∇· ( ρ v ) = ,∂ t ( ρ v ) + ∇· [ ρ vv + p ⋆ − BB ] = ρ Ω ( q Ω x ˆ x − ˆ z × v ) − ρ ∇ Φ ( z )ˆ z ∂ t B − ∇× ( v × B − η ∇× B ) = , ∇· B = , with a static gravitational potential Φ ( z ) = Ω z , the total pres-sure p ⋆ = p + B , and where all other symbols have their usualmeanings. The first source term in the momentum equation in-cludes the Coriolis acceleration and tidal force in the shearingbox approximation (see Gressel & Ziegler 2007; Stone & Gardiner2010, for implementation details). For the shear rate, q , we applythe Keplerian value of − / The simulations are semi-global in the sense that we only con-sider a small horizontal patch of a protoplanetary accretion discbut include the full vertical structure. Our fiducial stratified simu-lation covers ± H , to account for hydro-static stratification under the assumption of an isothermal equationof state. To accommodate su ffi ciently wide MRI-active regions inthe presence of a realistic dead zone, we use an increased box sizeof 5 . L z = H . Thisprovides an extra 3 to 4 H above and below the nominal dead zone,whose vertical extent either side of the midplane is ∼ H .Based on our previous work on permissible shearing box sizesfor studying the dynamics of embedded bodies in turbulent discs inpaper I, we chose a standard horizontal extent of 4 × H . Due tothe higher computational cost of the models including a dead zone,we adopt a somewhat smaller horizontal domain size of 3 × H ,yet still large enough to allow for the excitation of spiral densitywaves (Heinemann & Papaloizou 2009a,b). The box sizes and gridresolution for the di ff erent models are listed in Tab. 1.Because the ionisation model used to compute the resistivityinvolves chemical rate equations, we need to specify a unit sys-tem to convert between computational units and physical units forthe mass density and temperature. To allow for direct compari-son with other results, we scale our model such that it approxi-mates the widely used “minimum-mass” protosolar nebula model(Hayashi 1981). More specifically, we place our local box at a ra-dius r = h ≡ H / r = . Σ =
135 g cm − , a temperature T =
108 K,and an isothermal sound speed of c s =
667 m s − . Using our def-inition of H , the equilibrium vertical density profile is given by ρ ( z ) = ρ exp h − z / (2 H ) i , where ρ is the midplane density. domain [ H ] resolution XR SR CRA1 4 × , ± . × × − − − D1 3 × , ± . × × ◦ × ◦ D2 3 × , ± . × × × × ◦ B1 4 × × × × − − − Table 1.
Simulation parameters for the studied models. Columns 4-6 indi-cate the ionisation flux for X-rays (XR), short-lived radionuclides (SR), andcosmic rays (CR), respectively; open circles indicate nominal values, anda multiplication symbol followed by an integer gives the multiple of thenominal value adopted.
For the initial magnetic field, we apply the standard zero-net-fluxmagnetic configuration B z ( x ) ∼ sin(2 π x / L x ), with an initial plasmaparameter (ratio of thermal to magnetic pressure) β p ≃
50. To avoidlow β p in the far halo region, we attenuate B z ( z ) with height, as tokeep the plasma parameter roughly constant. At the same time, toobey the divergence constraint, we introduce a radial field B x ( x , z ),which is subsequently sheared out into an azimuthal field as thesystem evolves. The described configuration has proven to producea relatively smooth transition into developed turbulence, avoidingextreme field topologies in the linear growth phase.Following Turner & Sano (2008), we furthermore impose aweak additional B z net field, with an associated midplane β p = B x and B y near the top and bottom of the domain, which resultsin an accelerated wind in this region, clearing material and em-bedded field structures from the magnetically dominated halo re-gion. To compensate for the mass loss associated with the outflowboundaries and to enforce a stationary background, we continuallyre-instate the initial density profile in each grid cell by means of anartificial mass source term (cf. Hanasz et al. 2009).In the simulations adopting a dead zone due to insu ffi cientionisation we prescribe a locally and temporally varying mag-netic di ff usivity η . Some earlier simulations of dead zones byFleming & Stone (2003) and Oishi et al. (2007) prescribed a staticdi ff usivity profile derived from a physically motivated ionisationmodel. Here we further extend the realism of the di ff usivity modelby deriving a time-variable profile based on the gas column den- c (cid:13) , 1–18 Gressel, Nelson & Turner sity, the detailed procedure of which is described in the followingsection.After our disc models have achieved a steady turbulent state,we introduce 25 planetesimals into the disc midplanes. We treatthese planetesimals as non interacting test particles which evolveunder the full 3D influence of the gravitational field of the turbulentdisc. The disc gravity is softened on a length scale equal to the celldiagonal. ff usivity model Comparing di ff erent chemical reaction networks, Sano et al. (2000)and Ilgner & Nelson (2006a) found that the adsorption of free elec-trons and ions onto the surfaces of small dust grains plays an im-portant role in removing free charges from the gas phase. Theseauthors find that even tiny mass fractions of micron-sized grainscan significantly deplete charge carriers.Taking this dominance of small grains into account, wehere refrain from following the detailed non-equilibrium chem-istry in combination with a multi-fluid treatment as was done byTurner et al. (2007) and Ilgner & Nelson (2008). Instead we adopta simplified treatment of the gas-phase reactions. As a reason-able approximation to the more intricate modelling outlined above,we assume that recombination happens much faster than any dy-namical mixing timescale in our system, which seems warrantedin the presence of small grains (Sano et al. 2000; Ilgner & Nelson2006a,b). As a first step towards more realistic description, we al-low the magnetic di ff usivity, η = η ( x , t ), to vary spatially in allthree dimensions, and update it according to a look-up table derivedfrom the reaction network in model4 of Ilgner & Nelson (2006a).In particular, we assume dust grains of size 0 . µ m and with den-sity 3 g cm − and use a dust-to-gas mass ratio of 10 − . The assumedgas-phase Magnesium abundance is taken to be depleted by a factor10 − compared to its solar value. The input parameters to the tabu-lated di ff usivity are: (i) mass density, (ii) gas temperature, and (iii)the local ionisation rate ζ . All three quantities are evaluated on aper-grid-cell basis. To include the e ff ects of external irradiation on ζ , we compute column densities to both the upper and lower discsurfaces. There are several likely sources of ionising radiation in the vicinityof a newly born star. Turner & Drake (2009) have recently studiedhow contributions from stellar X-rays, radionuclides, and energeticprotons (interstellar cosmic rays, protons accelerated in the coronaeof the star and disc) influence the shape and extent of a possibledead zone. We have implemented their ionisation model, includingall ionisation sources. Because of uncertainties related to some ofthe proposed mechanisms, we take a conservative standpoint in thispaper and focus on stellar X-ray irradiation (XR), and interstellarcosmic rays (CRs) as the prime sources of ionisation.Observations of young stars in Orion show that T Tauri starX-ray luminosities vary by approximately four orders of magni-tude, with a median value of L XR ≃ × erg s − (Garmire et al.2000). Igea & Glassgold (1999) performed Monte-Carlo radiativetransfer calculations of the ionisation rate in a standard protostel-lar disc model, adopting the above value for L XR and assuming aplasma temperature of k T = ζ XR = . × − s − h e − Σ a / Σ XR + e − Σ b / Σ XR i r − , (1) Figure 1.
Vertical η profiles for runs D1 and D2, in terms of the magneticReynolds number, Rm = Ω H η ( z ) − , as given by our ionisation model. Thethick line corresponds to the adopted initial model (with Σ =
135 g cm − ).Thin lines on either side are for disc masses increased and reduced by fac-tors of two, and demonstrate the strong dependence on the mass column. Adead zone is expected for Rm ∼ <
10 (cf. Turner & Sano 2008). where r au is the position of our box in astronomical units, Σ a and Σ b , are the gas column densities above and below a given point,and Σ XR = . − is the absorption depth for the assumed spec-trum of X-rays. As detailed in Umebayashi & Nakano (2009), weprescribe the following vertical attenuation of interstellar CRs illu-minating the disc surfaces: ζ CR = × − s − e − Σ a / Σ CR + Σ a Σ CR ! − + . . . (2)Here, Σ CR =
96 g cm − is the cosmic ray attenuation depth esti-mated by Umebayashi & Nakano (1981), and the dots indicate thecorresponding contribution from the second column density Σ b .Because the numerically allowed di ff usive time step decreaseswith the grid spacing squared, adequately resolved dynamical sim-ulations can become prohibitive in the presence of high di ff usioncoe ffi cients. Given limited computational resources, we thereforedecided to confine the dynamical span in the magnetic di ff usivity η to a reasonable range. We primarily do this by including an am-bient ionisation due to the decay of short-lived radionuclides (SR),which we further enhance by a factor of ten compared to the fidu-cial value of ζ SR = . × − s − given in Turner & Drake (2009).As can be seen in Fig. 1, where we plot the resulting η profiles, thefloor in the magnetic Reynolds number due to the ambient level ofionisation is well below the expected threshold for the formationof a dead zone. This implies that we should not expect our resultsto be significantly a ff ected by enhancing the e ff ects of decayingradionuclides.We remark that, compared to previous work, our di ff usivityprofiles overall imply lower magnetic Reynolds numbers. Giventhat Fleming et al. (2000) have found a critical Rm c ≃ for sus-tained MRI in a zero-net-flux configuration, we expect that the am-plitude of an external net-flux will have a significant impact on thelevel of turbulence. Notwithstanding the moderate range in η , the numerical constraintgiven by the di ff usive propagation of information is one of the mainlimiting factors in our simulations. To avoid a degradation of the c (cid:13)000
96 g cm − is the cosmic ray attenuation depth esti-mated by Umebayashi & Nakano (1981), and the dots indicate thecorresponding contribution from the second column density Σ b .Because the numerically allowed di ff usive time step decreaseswith the grid spacing squared, adequately resolved dynamical sim-ulations can become prohibitive in the presence of high di ff usioncoe ffi cients. Given limited computational resources, we thereforedecided to confine the dynamical span in the magnetic di ff usivity η to a reasonable range. We primarily do this by including an am-bient ionisation due to the decay of short-lived radionuclides (SR),which we further enhance by a factor of ten compared to the fidu-cial value of ζ SR = . × − s − given in Turner & Drake (2009).As can be seen in Fig. 1, where we plot the resulting η profiles, thefloor in the magnetic Reynolds number due to the ambient level ofionisation is well below the expected threshold for the formationof a dead zone. This implies that we should not expect our resultsto be significantly a ff ected by enhancing the e ff ects of decayingradionuclides.We remark that, compared to previous work, our di ff usivityprofiles overall imply lower magnetic Reynolds numbers. Giventhat Fleming et al. (2000) have found a critical Rm c ≃ for sus-tained MRI in a zero-net-flux configuration, we expect that the am-plitude of an external net-flux will have a significant impact on thelevel of turbulence. Notwithstanding the moderate range in η , the numerical constraintgiven by the di ff usive propagation of information is one of the mainlimiting factors in our simulations. To avoid a degradation of the c (cid:13)000 , 1–18 ynamics of planetesimals in dead zones numerical accuracy in the MRI-active region by an excessively lownumerical time step in the dead zone, we choose to account for thetwo separate regimes in Rm. This is done by operator-splitting thedi ff usive term in the induction equation and applying a sub cycling-scheme for its update . By doing so, we are able to integrate thenon-di ff usive part of the MHD equations with a longer time step,enhancing the accuracy of the solution in regions of low magneticdi ff usivity. Properly resolving the growing modes of the magneto-rotationalinstability (MRI) with Godunov-type codes has been found to de-pend on the reconstruction strategy used and on the ability of theRiemann solver to capture the Alfv´enic mode (Balsara & Meyer2010). To improve the representation of discontinuities in our finitevolume scheme, we recently extended the nirvana - iii code with theHarten–Lax–van Leer Discontinuities (HLLD) approximate Rie-mann solver introduced by Miyoshi & Kusano (2005). To guaran-tee the required directional biasing of the electromotive force in-terpolation (cf. Flock et al. 2010), we have implemented and testedthe upwind reconstruction procedure of Gardiner & Stone (2005).To enhance the stability of our code in the strongly magne-tised corona, we gradually degrade the reconstruction order fromsecond- to first-order accuracy near the vertical domain boundaries,thus avoiding undershoots in the hydrodynamic state variables instrong shocks. ff usion To facilitate the study of a low-beta disc corona, Miller & Stone(2000) have introduced the concept of a so-called Alfv´en speed lim-iter, circumventing prohibitively high signal speeds in low-densityregions. Such a limiter can in principle be adopted for the approx-imate Riemann solvers we use. We chose, however, a di ff erent ap-proach and instead add an artificial mass di ff usion term to the equa-tions of mass and momentum conservation. The di ff usion coe ffi -cient is chosen to lie well below the truncation error in the bulkof the domain. In grid cells where the density contrast exceeds aspecified dynamic range, C dyn , the coe ffi cient is gradually adjustedaccording to D ( ρ ) = + C dyn ρρ ! − / (3)such that the grid Reynolds number approaches order unity – co-inciding with the stability limit for the explicit time integration ofthe di ff usive term. The transition function is chosen in a way thatthe inverse operation can be e ffi ciently implemented by means ofconsecutive square-root operations. We typically chose C dyn = ff usive fluxes are part of the finite vol-ume update, this approach does not violate mass conservation, andtherefore avoids the problems of enforcing an artificial floor valuein the density. In conclusion, our approach has proven to greatlybenefit the overall robustness without noticeably a ff ecting the solu-tion in the interior of the domain. Similarly to the limiters used byMiller & Stone (2000) and Johansen & Levin (2008), the Courant-Friedrichs-Lewy time-step constraint due to the fast magnetosonic We typically chose a sub-cycling ratio of 5 − Figure 2.
Visualisation of the turbulent flow structure. The colour codingindicates the vorticity, log |∇ × u | of the perturbed flow u = v − q Ω x ˆ y .While the turbulent regions above and below the dead zone show stronglyfolded vortex structures, the flow pattern near the midplane is dominated byshearing waves. The contrast in the vorticity amplitude is ∼ mode is significantly alleviated, allowing for economical use ofcomputing resources. In the following, we will present results from three di ff erent simula-tions (see Tab. 1 for details). To be able to make a quantitative com-parison with respect to the e ff ects of a dead zone, we performed afiducial stratified model with fully active MRI, referred to as modelA1. In general, this model agrees well with the unstratified sim-ulations presented in paper I. For our standard dead zone model,D1, we choose the reference ionisation rates of Turner & Drake(2009) along with a dust-to-gas mass ratio of 1 : 1000, i.e., account-ing for a modest depletion of the smallest grains by coagulation andsedimentation. The resulting di ff usivity profile is plotted in the leftpanel of Figure 1.The resulting dead zone in this model covers roughly ± H ,making it a reasonable proxy for what a realistic dead zone at5 au might look like. Gaseous nebulae around newly forming stars,however, are subject to substantial variations in X-ray irradiation(Garmire et al. 2000), so for our second dead zone model, D2, weincrease the stellar flux by a factor of twenty. We find the averagetotal thickness of the dead zone to be reduced by roughly 1 . One important result of the early simulations of layered accretiondiscs by Fleming & Stone (2003) was that the dead zone, despite itsname, retains a non-negligible level of Reynolds stresses, namely in c (cid:13) , 1–18 Gressel, Nelson & Turner orbits h α i σ Γ [ cm s − ] τ c [2 π Ω − ] C σ ( ∆ x ) [ H ] C σ ( e ) [ H / r ]A1 217 0.0105 0.45 × × − × − D1 224 0.0038 0.06 × × − × − D2 223 0.0051 0.13 × × − × − B1 ∗
505 0 .
05 0.74 × × − × − ∗ torque corrected for 2D /
3D evaluation (cf. Fig. 8 in Nelson & Gressel 2010).
Table 2.
Overview of simulation results. Gravita-tional torques Γ represent the standard deviation σ of a normal distribution (see Fig. 8), Coherence times τ c are from a fit to the torque ACF (see Fig. 9). Therandom walk amplitudes C σ ( ∆ x ) and C σ ( e ) refer tothe dispersion in the displacement and eccentricity ofa swarm of particles. Figure 3.
Comparison between the box averaged turbulent Maxwell- andReynolds stresses. Values are normalised by the gas pressure in the discmidplane, p . The time-averaged sums of these, i.e., the α parameters, arelisted in column three of Tab. 2. the form of waves that are excited in the active layers. These resid-ual motions can clearly be seen in Fig. 2, where we have visualisedthe flow structure in terms of the logarithmic vorticity. The colourcoding exhibits very distinct patterns in the two regions: while theMRI active layers show folded vortex features characteristic of de-veloped turbulence, the dead zone region is clearly non-turbulent,and dominated by sheared density waves with k R ≫ k φ . At theinterface between the two zones, one can furthermore see struc-tures indicating some level of turbulent mixing into the dead zone.Because our model assumes that recombination on the surface ofsmall dust grains is e ffi cient (leading to a recombination time scalemuch smaller than the dynamical mixing time scale) this does not,however, have an e ff ect on the extent of the dead zone.Our primary aim in this paper is to examine statistically thedynamical evolution of planetesimals embedded in turbulent discswith and without dead zones. In order to do this we require a quasi-stationary model of turbulence. As can be seen in Fig. 3, where weplot the time variation of the box-averaged transport coe ffi cients,our models fulfil this requirement of a stationary mean with (ad-mittedly strong) superposed fluctuations. Because the bulk of thetransport is associated with the MRI-active regions, the turbulentstresses only seem to depend weakly on the actual width of thedead zone. We remark that due to the di ff erent field configuration,the overall strength of the turbulence is somewhat weaker than inour previous non-stratified models presented in paper I. This is re-lated to buoyant loss of magnetic flux from the disc in stratifiedsimulations, and is consistent with recent results in the literature(cf. Davis et al. 2010, and references therein). Figure 4 shows space-time plots of the mean toroidal magnetic field¯ B y ( z , t ), as well as the total Maxwell stress for the fully active modelA1. Probably the most striking features in this plot are the peri- odic “butterfly” patterns in the large-scale magnetic field, whichare commonly observed in this type of simulation (e.g. Davis et al.2010; Flaig et al. 2010). Because the upwelling field structures arenot associated with a bulk motion of the flow, the likely explana-tion for the emerging patterns is a mean-field dynamo as recentlyre-investigated by Gressel (2010). Even though we see quite strongfluctuations in the total Maxwell stress, the overall state is quasistationary. Note that the two most prominent peaks of the totalMaxwell stress (i.e. around 140-145 and 155-160 orbits) are (a)due to a strong mean field, and (b) confined to one half of the box.The situation of magnetic activity being dominant in one half ofthe disc has been observed before (see e.g. fig. 7 in Miller & Stone2000), and might be explained by equally strong quadrupolar anddipolar contributions. These cancel on one side and, at the sametime, add-up on the opposite side, resulting in the observed lop-sided appearance. In fact, we cannot define a clear parity of thefield in Fig. 4, which implies that the dipolar and quadrupolar dy-namo modes seem to possess similar growth rates.In the following discussion, we will focus on the vertical struc-ture of the disc and how it is shaped by the presence of a dead zone.For this we will show the time-averaged vertical profiles of vari-ous quantities. To demonstrate that this is a meaningful procedure,we provide space-time diagrams of three representative quantities(Fig. 5) for model D1. Again, the mean toroidal field generally isof irregular parity, i.e., neither quadrupolar nor dipolar symmetry isseen to prevail. Whether the concept of a global parity is relevantin the presence of a highly di ff usive “insulating” layer remains amatter of discussion. It is interesting to note, however, that MRIchannel modes – which are a non-linear solution to local net-fluxMRI simulations – represent global modes even when stratificationis included (Latter et al. 2010). In fact, we see some indication ofchannels in our simulations, most prominently in the density (com-pare the lower panel in Fig. 5 with Latter et al. 2010, fig. 3).Quite remarkably, the number of field reversals related to thedynamo-cycle seems to be largely una ff ected by the presence of thedead zone. As mentioned above, one main motivation in explainingthe rising patterns by a mean-field dynamo, was the absence of abulk motion near the midplane – implying that the pattern speedcannot be explained by the Parker instability (Shi et al. 2010). Thekey in explaining the upward motion was a negative buoyant α ef-fect near the midplane (R¨udiger & Pipin 2000; Brandenburg 2008).With the pattern in the disc halo being unchanged in the pres-ence of a dead zone, this picture possibly has to be amended.Naively, for a positive α e ff ect in the halo (as found by Gressel2010), one would not expect an outward butterfly diagram fromclassical theory – demonstrating the need for non-linear feedbackboth in the diagnostics (Rheinhardt & Brandenburg 2010) as wellas in the modelling (Brandenburg, Candelaresi & Chatterjee 2009;Courvoisier, Hughes & Proctor 2010).We note however, that because of the presence of the weak net If not stated otherwise, we average over the interval [20,210] orbits.c (cid:13) , 1–18 ynamics of planetesimals in dead zones Figure 4.
Spatio-temporal evolution of (a) mean toroidal field, (b) total Maxwell stress, for the fully active model A1.
Figure 5.
Same as Fig. 4, but for the dead zone model D1, and including an extra panel showing the relative density fluctuations (logarithmic). The horizontalfield exhibits the same dynamo cycles in the MRI active layers. As already seen by Turner & Sano (2008), the toroidal field leaks into the di ff usively dominatedmidplane region, contributing to the stresses. Despite the absence of strong turbulent fields, density fluctuations reach a level of 10-20% within the dead zone. vertical flux in our current simulations, they are not strictly compa-rable to our previous work. Due to the net flux, a possible explana-tion of the unchanged migration direction might be the presence ofprevailing channel modes. These were shown to produce a strongnegative helicity (cf. fig. 4 in Gressel 2010), compatible with theupward motion of the wave patterns .Now focusing our attention to the actual dead zone region, wesee that a significant amount of azimuthal (and in fact, radial) fielddi ff uses in to the low-Rm region around the midplane. This leak-age has first been seen in simulations by Turner et al. (2007) and This was already suggested as an alternative scenario in the mentionedpaper.
Turner & Sano (2008), which are very similar to ours and also con-tain a weak net vertical field. In contrast to this, the simulationsby Fleming & Stone (2003), Oishi, Mac Low & Menou (2007) andIlgner & Nelson (2008), which apply a zero-net-flux configura-tion, do not show such an e ff ect. It seems peculiar, however, thatthe presence of a net flux should play a role in this context. AsIlgner & Nelson (2008) illustrate in their fig. 2, the assumed di ff u-sivity profiles of Fleming & Stone (2003) result in moderately highmagnetic Reynolds numbers near the midplane. On the contrary,in the model of Turner, Sano & Dziourkevitch (2007) and also inour model (cf. Fig 1) the magnetic Reynolds number reaches muchlower values of Rm <
10 for | z | < H , leading to a much reduced c (cid:13) , 1–18 Gressel, Nelson & Turner
Figure 6.
Maxwell stresses for the di ff erent models. Thick lines indicatethe total stress ∝ − h B R B φ i , thin lines represent the turbulent contribution ∝ − h B ′ R B ′ φ i , where B ′ = B − h B i . The vertically integrated mass accretionrates (see rhs axis) for the models D2, D1 and A1, are 9 . × − , 7 . × − , and 1 . × − M ⊙ yr − , respectively. Note the unexpected trend inthe depth between models D1 and D2, as opposed to Fig. 7. di ff usion time scale, possibly facilitating the leakage of large-scalefield from the disc halo into the dead zone.Finally, we note that relative density fluctuations remain at alevel of up to ten per cent within the dead zone – as seen in the low-ermost panel of Fig. 5. Episodes of lower perturbance (as indicatedby lighter colours) are clearly correlated with strong azimuthalfield, and notably with a negative Maxwell stress at the interfacebetween the dead and active zones. As already mentioned above,we see indications of persistent channel modes (dark streaks above | z | ∼ > H ), which are not unexpected given the magnetic configura-tion. It might be argued that such features are overly pronouncedat the current numerical resolution, and might be destroyed byparasitic instabilities (Pessah & Goodman 2009; Latter et al. 2009,2010) when the resolution is further increased. Figure 6 shows time-averaged vertical profiles of the R φ compo-nents of the Maxwell stress tensor. Because we use an isothermalequation of state, resulting in a constant sound speed, c s , these pro-files translate directly into a fractional mass accretion rate as in-dicated by the right-hand axis. The values on this axis are givenby the expression H d ˙ M / d z = π c s H ρ ( z ) T r φ / p ( z ), which can bederived from the usual approximation for the mass accretion ratein a steady disc ˙ M = πν Σ , but relaxing the assumption that thee ff ective viscous stresses are independent of z .The vertically integrated mass flow rates generated by themagnetic stresses are 9 . × − and 7 . × − M ⊙ yr − for runsD2 and D1, respectively. This similarity is expected because theheights of the dead and active zones about the midplane in eachmodel are similar ( ∼ . H and 2 H , respectively, for the dead zonesin models D2 and D1). The mass flow rate in the fully active run A1is 1 . × − M ⊙ yr − , nearly twice that observed in the model withthe largest dead zone. Given that the bulk of the mass in the disclies within the dead zone, where the Maxwell stresses for modelsD2 and D1 are factors of ∼ − and 10 − smaller than for modelA1, the bulk of the transport in models D1 and D2 happens in theactive layers because the stresses (and e ff ective viscosity coe ffi cient α ) remain large there.Figure 6 shows that the Maxwell stresses in the midplaneregions are dominated by large-scale fields in runs D1 and D2 Figure 7.
Same as Fig. 6, but for the Reynolds stress h ρ v R δ v φ i . Thevertically-integrated mass accretion rates are 2 . × − , 1 . × − , and3 . × − M ⊙ yr − for models D2, D1, and A1, respectively. – even more so than in the halo regions of the fully active runA1. This result agrees broadly with recent simulations (cf. fig. 4in Turner & Sano 2008) which included similar physical e ff ects.While the turbulent contribution to the Maxwell stress (thin lines)falls-o ff to ≃ − in the dead zone for both models D1 and D2,one can clearly see that the e ff ect of the greater dead zone widthin model D1 is compensated for by a shallower valley in the totalMaxwell stress (thick lines). The inverse top-hat shape of the pro-files is likely related to large-scale field di ff using into regions ofintermediate Rm, as seen in the topmost panel of Figure 5. Giventhat the transition region is closer to the low- β halo for model D1(with a wider dead zone), it seems plausible that the strength of theregular field is larger in this case. Clearly, the contribution due tothe fluctuating fields becomes negligible in both cases.Considering that, for model D1, the floor value of the Maxwellstress is of the same magnitude as the Reynolds stress (cf. Fig. 7),such large-scale stresses might play an important role in mediat-ing mass transport within the dead zone – making a more detailedfuture study of the discussed trend worthwhile, particularly in thecontext of global disc models, where it may be possible to observelarge scale accretion flows in dead zones.In Figure 7, we plot vertical profiles of the time-averagedReynolds stress. As expected, the lines agree very well in the MRI-active regions. Comparing the models D1 and D2, we now see acorrespondence between the dead zone width and depth, as ex-pected. It appears that larger and more massive active zones lyingabove and below the dead zone are more able to excite higher am-plitude sound waves which are able to propagate into the dead zone,generating a correspondingly larger Reynolds stress there.Based on two runs with di ff erent (zero net flux) field topolo-gies, Fleming & Stone (2003) concluded that “the Reynolds stressin the midplane does not seem to depend on the size of the deadzone but rather the amplitude of the turbulence in the active layers”.This is clearly not the case for our models D1 and D2, which showquite di ff erent Reynolds stresses in the dead zone, while havingalmost identical levels of Maxwell stresses – both in the total andfluctuating part – within the active region (see Figs. 6 and 7, respec-tively). While we reckon that the larger Reynolds stress is relatedto the larger mass of the active layer in model D2, this discrepancymay also be due to the presence of large-scale fields di ff using intothe midplane region in our simulations. Note that Fleming & Stone(2003) did not mention such fields in their discussion. c (cid:13) , 1–18 ynamics of planetesimals in dead zones Figure 8.
Torque histograms for models A1, D1 and D2. The grey lines indicate fitted normal distributions, for which standard deviations are given (note thedi ff erent scale on the abscissae). While the torques within the dead zone are well approximated by a Gaussian, the distribution in model A1 develops somedeviations for large torques. As first suggested by the global simulations of Nelson (2005), thestochastic gravitational forces associated with density fluctuationsfrom developed MRI turbulence have the potential to limit severelythe growth of planetesimals. In paper I, we have confirmed thisoriginal finding in the framework of global simulations and localsimulations with large enough box sizes, but in which the verticalcomponent of gravity was neglected. In the following section, weextend this work to the case of stratified discs with and withouta dead zone. To facilitate a direct comparison, we here focus onthe case of a moderately strong net vertical field. In this respect,the results should be interpreted as upper limits. To obtain lowerlimits for the torque amplitude, models with nominal dead zonesand applying a zero-net-flux configuration will be required.
Figure 8 shows distribution functions of gravitational torques alongthe y -direction in our simulation box. The histograms are computedfrom the time series of a set of 25 particles in each simulation run.Unlike Oishi et al. (2007), we do not observe transients in the timeseries, and the torques are indeed quasi-stationary in the interval[20,210] orbits, which has been used to obtain each histogram. Asrequired for a simplified treatment in terms of a stochastic MonteCarlo model, the histograms are well represented by a normal dis-tribution. Notably, there are moderate power-law tails in the fullyactive run A1, indicating some level of non-Gaussian fluctuations(also cf. fig. 7 in paper I). The amplitude of the torques in thedi ff erent runs roughly scales with the square-root of the midplaneReynolds stress. This arises because the density fluctuations scalelinearly with the velocity fluctuations, and coincides with the scal-ing found in section 3.3 of Baruteau & Lin (2010).Apart from the torque amplitude, the ability to stochasticallyamplify particle motions depends on the degree of coherence inthe fluid patterns. Building on our work from paper I, we here ap-ply the same formalism and study temporal torque autocorrelationfunctions (ACFs). We use the fitting formula introduced in section3.3 of paper I, i.e. S Γ ( τ ) = [(1 − a ) + a cos(2 π ω τ )] e − τ/τ c , (4)with free parameters, a , ω , and τ c . As discussed in paper I, we as-sume two components to better fit the shape of the ACF: (i) an ex-ponential decay representing the truly stochastic part of the torquetime series; (ii) a modulation due to “wavelike” behaviour, i.e., a Figure 9.
Torque autocorrelation function (ACF) for models A1 and D1.The black lines indicate model fits (see text), showing that motions are co-herent over approximately one third of the orbital time scale. The ACF formodel D2 is very similar to D1 and is hence not shown. negative dip in the ACF reflecting the fact that density enhance-ments swashing over the region of influence will produce torquesof opposite signs when approaching and receding the test particle.Such a feature is also seen in the corresponding Figs. 9 and 15 ofOishi et al. (2007) and Yang et al. (2009), respectively, complicat-ing the interpretation of the results.In paper I, we have discovered certain trends in the ampli-tude of the latter e ff ect, which are probably related to the finite sizeof the computational domain. Because the periodic box in somerespect acts as a resonator, the wavelike component probably ap-pears too pronounced in local simulations. This demonstrates thatone does well keeping in mind the limitations of local boxes interms of the allowed dynamics (also see Regev & Umurhan 2008).The parameter of interest, of course, is the correlation time τ c ofthe stochastic perturbation. While waves represent ordered motion,it is this stochastic part that ultimately produces the random-walkamplification of the particle motions.In the following, we rely on the approach taken in Eqn. (4)providing an e ff ective deconvolution of the two e ff ects. The result-ing fits for models A1 and D1 are shown in Fig. 9. Unlike modelB1 (cf. fig. 9 in paper I), model A1 shows a moderate level of mod-ulation. Given the horizontal box dimensions of 4 × H , this issomewhat unexpected as we observed a consistent trend towardsweaker modulation for larger box sizes in paper I. This being said,the fitted coherence time of τ c = .
30 in model A1 is in strikingagreement with our earlier unstratified model. The torque ACF of c (cid:13) , 1–18 Gressel, Nelson & Turner
Figure 10.
Comparison ofthree di ff erent model ACFs asgiven by Eqn. (4). All curveshave the same coherence time τ c = . ff erent zero-crossings. the dead zone model D1 shows a comparable amplitude in the mod-ulation and has an almost identical value of τ c = . τ c = . ff erentmodulation of the exponential decay as exemplified in Figure 10.All curves have a common τ c = .
3. Moreover, the first two curveshave a 60% wavelike modulations with frequency ω as indicated,while the third is a simple exponential. Owing to the confusion withrespect to the definition of the coherence time, it seems worthwhileto further study how the spectrum of modes a ff ects the apparentcoherence. Nevertheless, given the actual type of forcing used inBaruteau & Lin (2010), i.e. a spectrum of wavelike motions, onemay be surprised how similar the ACFs indeed look. Going back toFig. 2, it however becomes obvious that it is not the turbulence weneed to approximate but only the spiral density waves with uncor-related phases that it creates.The corresponding ACF for model D2 is omitted here, as itis very similar to D1. The fitted coherence time for model D2 is τ c = .
27, in good agreement with the other two runs. We thereforeconclude that the presence (and depth) of a dead zone does havelittle influence on the temporal torque statistics and merely a ff ectsthe amplitude of the stirring process. To study the excitation of the velocity dispersion and the radialdi ff usion of embedded boulders, planetesimals, and protoplanets,in each run we integrated the trajectories of a swarm of 25 non-interacting test particles subject to perturbations from the gravita-tional potential of the gas disc. Unlike in paper I, we here focuson larger bodies and neglect the e ff ects of aerodynamic interaction.This is warranted for solids with radii above ≃
100 m (cf. paper Iwhere the results for planetesimals with sizes in the range 100 m -10 km were found to be similar for run times on the order of a fewhundred planetesimal orbits). The upper limit for the size rangeconsidered here is given by the constraint that perturbations of thedisc remain weak such that spiral wave excitation and gap-openingcan be ignored. This is the case for planetesimals and small proto-planets, where feedback onto the disc can be ignored.Although we do not consider in detail the dynamics of metre-sized boulders that are strongly coupled to the gas via drag forces inthis paper, we recall from paper I that such bodies quickly achieverandom velocities very similar to the turbulent velocity field of thegas. The midplane r.m.s. turbulent velocity for the fully turbulentmodel A1 is found to be v turb = (76 . ± .
5) m s − , and for modelD1 v turb = (17 . ± .
6) m s − . For model D2 v turb = (28 . ± .
6) m s − . Figure 11.
Random-walk eccentricity growth for the dispersion, σ ( e ), ofsets of 25 particles in each run. The expected √ t behaviour is indicatedby the dashed line. The rhs axis gives the corresponding radial velocitydispersion, leading to disruptive relative velocities for the fully active runA1. The amplitudes for the dead zone runs are lower by a factor of 10-20. We now consider the dynamics of larger planetesimals. Theexcitation of the eccentricity (and equivalently the radial velocitydispersion) can be seen in Fig. 11, where we plot the r.m.s. val-ues of these quantities averaged over the planetesimal swarms as afunction of time for the runs A1, D1 and D2. The time history isconsistent with a random-walk behaviour, indicated by the dashedline representing an e ( t ) ∼ √ t dependence. Expressing the timeevolution of the eccentricity according to e ( t ) / h = C σ ( e ) √ t , (5)where h = H / r , and the time is measured in local orbits, we obtainthe following values for the amplitude factors, C σ ( e ) (which arealso listed in Tab. 2): 2 . × − , 2 . × − , and 1 . × − formodels A1, D2, and D1 respectively. Despite the di ff erent levelsof turbulence (as reflected in both α and Γ ), the values in Tab. 2show that the stirring amplitude of model A1 agrees quite well withthe unstratified model from paper I. The amplitudes from the deadzone runs D2 and D1 are lower by a factor of ≃
11 and ≃ ff ect on the strength of eccentricity stirring. As we have only been able to run our simulations for about 200orbits, it is important to consider the long term evolution by exam-ining the magnitude of the equilibrium eccentricity that is expectedto arise through the balance between turbulent eccentricity excita-tion and gas drag and / or collisional damping (Ida et al. 2008). Thequestion of whether or not collisional growth of planetesimals mayoccur in the di ff erent disc models is determined by the velocity dis-persions obtained relative to the disruption / erosion velocity thresh-olds, and this point is addressed in the discussion that follows. We have not included gas drag or collisional damping in the sim-ulations presented in this paper, but gas drag was included in theruns presented in paper I, and they showed that over run times ofa few hundred orbits the gas drag has little influence on the eccen-tricity for planetesimals with radii R p >
100 m. We adopt the usualexpression for the gas drag force (Weidenschilling 1977) c (cid:13)000
100 m. We adopt the usualexpression for the gas drag force (Weidenschilling 1977) c (cid:13)000 , 1–18 ynamics of planetesimals in dead zones Model size v disp [ m s − ] v disp [ m s − ] v crit [ m s − ] v crit [ m s − ] v crit [ m s − ] t eq [ yr ](gas drag) (collisional) (weak) (strong) (rubble pile)A1 100 m 73.8 79.9 1.0 7.7 - 1 . × . ×
10 km 342.8 795.0 26.3 137.0 58.5 4 . ×
100 km 738.2 2514.3 262.6 1368.0 585.4 1 . × D1 100 m 11.0 4.5 1.0 7.7 - 2 . × . ×
10 km 51.0 45.4 26.3 137.0 58.5 2 . ×
100 km 110.0 143.5 262.6 1368.0 585.4 10 . × D2 100 m 15.2 7.4 1.0 7.7 - 2 . × . ×
10 km 70.5 73.9 26.3 137.0 58.5 1 . ×
100 km 153.4 233.6 262.6 1368.0 585.4 7 . × Table 3.
Results for the equilibrium velocity dispersions obtained for each model as a function of the planetesimal sizes, and as a function of the dampingmechanism. Also tabulated is the time required to grow to the smaller of the equilibrium velocity dispersion values given for each model each size. F drag = C D π R ρ | v p − v g | ( v p − v g ) (6)where C D is the drag coe ffi cient. For making simple estimateswe take the value C D = .
44, appropriate for larger planetesi-mals for which the Reynolds number of the flow around the body R e > t , is measured in local orbits. Following the same proce-dure adopted in paper I, and calculating the equilibrium velocitydispersion by equating the eccentricity growth time scale and thegas drag-induced damping time scale, we obtain v disp = C σ ( e ) h R p ρ p v C D ρ ! / (7)where R p and ρ p are the planetesimal radius and density, and v k isthe Keplerian velocity. Note that the form taken by (7) now assumesthat the system evolution time given in (5) is expressed in seconds.The equilibrium velocity dispersions obtained using (7) for eachdisc model and planetesimal size are listed in the third column ofTab. 3. The e ff ects of collisional damping, due to inelastic collisions be-tween planetesimals, may be estimated by noting that the colli-sional damping time scale is given by the product of the collisiontime scale and the number of collisions required to damp out thekinetic energy of random motion for a typical planetesimal. Thisapproach is similar to that adopted by Ida et al. (2008). The colli-sion time (neglecting the e ff ect of gravitational focusing) is givenby t coll = n p π R v disp , (8)where n p is the number density of planetesimals. This may be ap-proximated by n p = Σ p / (2 H p m p ), where H p = v disp / Ω k and Σ p isthe mass surface density of planetesimals. The coe ffi cient of resti-tution is given by C R = v f / v i , where v i and v f are the initial and finalvelocities associated with a collision. The number of collisions re-quired to damp v disp is therefore N coll = − C R . (9) Combining Equations (8) and (9) we obtain the collision dampingtime: τ c − damp = R p ρ p Σ p Ω k − C R ! . (10)Equating the collisional damping time with the eccentricity growthtime associated with (5) – and given explicitly by equation (16) inpaper I – yields an expression for the equilibrium velocity disper-sion v disp = s R p ρ p C σ ( e ) h v k Σ p Ω k − C R ! . (11)The value to be adopted for the coe ffi cient of restitution, C R ,depends on the material composition, the impact velocity, andwhether or not the planetesimals may be considered to be mono-lithic bodies or rubble piles held together by self-gravity alone.Adopting the velocity dependent formula from Bridges et al.(1984) C R = MIN ( . (cid:18) v disp − (cid:19) − . , ) , (12)we find that C R ≃ v disp , that we obtain in the simulations. We therefore adopt C R = v disp . The equilibrium velocity dis-persions obtained using (11) for each disc model and planetesimalsize are listed in the fourth column of Table 3. In obtaining thesevalues we assume that Σ p = × ( Σ / /
250 lower than the gas surface density, but isaugmented by a factor of 4 beyond the ice-line. We further assumethat all disc solids are incorporated within planetesimals of size R p for each size that we consider. The consequences of the equilibrium velocity dispersions obtainedfrom equations (7) and (11) for planetary growth can only be deter-mined by comparing them with the disruption / erosion thresholdsfor colliding bodies (Benz & Asphaug 1999; Stewart & Leinhardt2009), and with the escape velocities associated with the planetes-imals. In a recent study, Stewart & Leinhardt (2009) present a uni-versal law for collision outcomes in the form c (cid:13) , 1–18 Gressel, Nelson & Turner M lr M tot = − Q R Q D , (13)where M lr is the mass of the largest post-collision remnant, M tot is the total mass of the colliding objects M + M , and Q R is thereduced mass kinetic energy normalised by the total mass Q R = M M M V . (14)Here V I is the impact velocity. For accretion to occur during a col-lision between equal sized bodies, we require M lr / M tot > /
2, orequivalently Q R / Q D <
1, such that Q D is the collisional disruptionor erosion threshold.The value of Q D is sensitive to factors that influence the energyand momentum coupling between colliding bodies (e.g. impact ve-locity, strength and porosity). Stewart & Leinhardt (2009) fit resultsfrom their numerical simulations and data in the literature using theexpression Q D = h q s R µ m / (3 − φ m )12 + q G R µ m i V (2 − µ m )I , (15)where q s , q G , φ m and µ m are parameters. R is the spherical ra-dius of the combined mass, assuming ρ p = − . Since we use ρ p = − , and consider collisions between equal-sized bodiesonly, we have R = / R p . The first term on the right hand sideof (15) represents the strength regime, while the second term rep-resents the gravity regime. Small bodies are supported by materialstrength, which decreases as the planetesimal size increases. Onthe other hand, gravity increases in importance with growing plan-etesimal size – with the transition between the strength and gravityregimes occurring at R p ≈
100 m. Stewart & Leinhardt (2009) de-rive the following values for the above parameters for weak aggre-gates (weak rock): q s = q G = − , µ m = . φ m =
7. Forstrong rocks they use basalt laboratory data and modelling resultsfrom Benz & Asphaug (1999) to obtain: q s = × , q G = − , µ m = . φ m =
8. In the limit of large planetesimals, the disruptioncurve can be best fit by their results for colliding rubble piles, forwhich Q D = . × − R . Given the uncertainties associated withthe structure and material strength of planetesimals, we tabulate thedisruption velocities, v crit , for both weak and strong planetesimalsin the fifth and sixth columns of Tab. 3. For R p >
10 km, we tabulatethe disruption velocities for rubble piles in the seventh column.
We now consider the results for the fully turbulent model A1. Theequilibrium velocity dispersions corresponding to gas drag damp-ing are listed in the third column of Tab. 3 for 100 m, 1 km, 10 kmand 100 km-sized planetesimals (assuming a density ρ p = − ). These values of v disp are smaller than those presented in pa-per I by approximately a factor of two, because in that paper weadopted a slightly more massive disc model, a larger planetesimaldensity ρ p = − , and we utilised cylindrical disc models thatgive rise to enhanced stirring of planetesimals relative to the full3D simulations considered here.Comparing the values of v disp in the third and fourth columns,we see that for all planetesimal sizes damping due to gas drag dom-inates over collisional damping, so it is the values in the third col-umn that are closest to the true equilibrium values obtained whenall sources of damping act simultaneously. We see that v disp > v crit for planetesimals composed of both ‘weak’ and ‘strong’ rock for R p
10 km, and for R p =
100 km we see that v disp exceeds thedisruption / erosion threshold for rubble piles. We conclude that if planetesimals reach their equilibrium values for v disp , then mutualcollisions will lead to their destruction.Estimated times for v disp to grow to the equilibrium values viaa random-walk ( v disp ∼ √ t ) are listed in the eighth column of Tab. 3.The most favourable models for the incremental formation of plan-etesimals at 5 au in a laminar disc suggest formation times of afew × yr (Weidenschilling 2000). Formation times in a turbu-lent disc exceed this because the vertical settling of solids is re-duced (Brauer et al. 2008). The random velocity growth times inTab. 3 thus indicate that if planetesimals were able to form incre-mentally in the disc that we simulate, then they would always havea velocity dispersion equal to the equilibrium value. But, the timeto reach the disruption velocity v crit = . − for strong 1 kmaggregates is only 800 yr, much shorter than the formation time. Itis clear that planetesimals cannot form and grow by a process ofcollisional agglomeration in a fully turbulent disc similar to modelA1. We first discuss model D1, which has the deeper dead zone of thetwo dead zone simulations. The values of v disp listed in Tab. 3 showthat for planetesimals with R p
10 km collisional damping domi-nates gas drag in setting the equilibrium velocity dispersion. Onlyfor R p =
100 km is gas drag more important. The transition fromgas drag dominated damping in model A1 to collisional damping inmodel D1 occurs because of the di ff erent functional dependencieson the eccentricity excitation coe ffi cient, C σ ( e ) in equations (7) and(11).The equilibrium v disp for all planetesimal sizes lie between thedisruption / erosion thresholds, v crit , for weak and strong aggregates.For the larger planetesimals with R p =
10 and 100 km we see that v disp < v crit for rubble piles. These results suggest that collisionsbetween planetesimals in a dead zone can lead to growth rather thandestruction, with the outcome depending on the material strength.The times required for v disp to reach the equilibrium values arelisted in the eighth column of Tab. 3. Given that the time requiredfor incremental formation of km-sized planetesimals in a disc withmodest turbulence is likely to exceed 10 yr at 5 au, the time scalesfor the growth of v disp suggest that planetesimals growing througha process of particle sticking will always have velocity dispersionsclose to the equilibrium values. The fact that these values are belowthe disruption thresholds for strong aggregates suggests that colli-sional growth may still be possible, albeit at a slower rate than in alaminar disc.If km-sized planetesimals can form, further accretion isnormally expected to arise via runaway growth, leading tothe formation of ∼ ∼ yr(Wetherill & Stewart 1993; Kenyon & Bromley 2009). Runawaygrowth requires the velocity dispersion to be lower than the es-cape velocity from the accreting bodies, and for an internal density ρ p = − this is given by v esc = (cid:16) R p (cid:17) m s − . A typical10 km planetesimal in model D1 forming via incremental growthwill be excited to v disp >
10 m s − during its formation, so turbu-lent stirring will prevent runaway growth from occurring for bodiesof this size. Because v disp remains below the disruption thresholdfor rubble piles, however, collisions can still lead to growth, butat a substantially slower rate than occurs during runaway growth.We discuss the implications of our results for the more rapid plan-etesimal formation scenarios presented by Cuzzi et al. (2008) andJohansen et al. (2007) in Section 7.Similar conclusions may be drawn from the results of model c (cid:13)000
10 m s − during its formation, so turbu-lent stirring will prevent runaway growth from occurring for bodiesof this size. Because v disp remains below the disruption thresholdfor rubble piles, however, collisions can still lead to growth, butat a substantially slower rate than occurs during runaway growth.We discuss the implications of our results for the more rapid plan-etesimal formation scenarios presented by Cuzzi et al. (2008) andJohansen et al. (2007) in Section 7.Similar conclusions may be drawn from the results of model c (cid:13)000 , 1–18 ynamics of planetesimals in dead zones Figure 12.
Dispersion of the radial displacement ∆ x of a swarm of 25 testparticles for run A1. The growing spread in the particles’ separation is wellapproximated by a random-walk as illustrated by the over-plotted curve. D2 as we have drawn for model D1. But, having a deeper deadzone, as in model D1, is clearly preferable for planetesimal forma-tion due to the fact we find that the velocity dispersion is smaller inthat case.
It was shown in paper I that the radial drift of planetesimals of size R p >
10 km over typical disc life times due to gas drag is small. In-stead, the di ff usion of planetesimal semimajor axes is dominated bythe stochastic gravitational torques exerted on the planetesimals bythe turbulent disc. This implies that material of common origin andcomposition (e.g. planetesimals which form at a particular radiallocation in the disc), in time will spread out over a certain region insemimajor axis. The degree of spreading is determined by the disclife time and the strength of the stochastic torques.In paper I, we examined the degree of spreading induced bythe stochastic torques in a fully turbulent disc model similar inmass to the minimum mass solar nebula model, and showed thatover a putative disc life time of 5 Myr, planetesimals embedded inthe current asteroid-belt region would spread inward and outwardover a distance of ∼ . ff erent taxonomic classes of asteroidsare distributed as a function of heliocentric distance has becomemore complex in recent years (Moth´e-Diniz et al. 2003), the earlierwork of Gradie & Tedesco (1982) shows a clear trend of increas-ing volatile content as a function of heliocentric distance. In thislatter study, the distribution of taxonomic types is consistent witha picture of in situ formation, followed by radial di ff usion over adistance of ≃ . Figure 13.
Model comparison for the random-walk behaviour shown inFig. 12 above, but plotted logarithmically. The second axis gives the ab-solute dispersion ∆ a at r = In Sect. 4, we have demonstrated that the gravitational torquesacting on particles can be represented by a normal distribution, andtheir temporal correlations possess a finite coherence time. Thesetwo properties allow us to estimate particle di ff usion based on aFokker-Planck equation. As discussed in detail in section § j . The timescale for di ff usion of particle angular momenta due to stochastictorques is t j = ( ∆ j ) D j (16)where D j is the di ff usion coe ffi cient and ∆ j is the change in j . Thedi ff usion coe ffi cient D j may be approximated by D j = σ Γ τ c , where σ Γ is the r.m.s. of the fluctuating specific torques, and τ c is thecorrelation time associated with the fluctuating torques.The change in specific angular momentum obtained after anevolution time of t is given by ∆ j = p D j t . (17)Noting that small changes in specific angular momentum are re-lated to small changes in the semi-major axis, a , by the expression ∆ jj = ∆ a a (18)the change in semimajor axis over an evolution time t is: ∆ a = a p D j tj . (19)We can now examine the degree of agreement between thelevel of particle di ff usion observed directly in the simulations, andpresented in Fig. 13, and predictions based on (19). The radial lo-cation of the shearing box in our simulations is assumed to be 5 au,and simulation run times are t ≃
200 local orbits. The value of D j for each model may be computed using the expression D j = σ Γ τ c ,and values for σ Γ , expressed in c.g.s. units, may be read o ff Fig. 8for each model. The corresponding estimates of τ c are listed inTab. 2 (along with the r.m.s. specific torques, σ Γ ).We note that the time evolution of the r.m.s. radial displace-ment obtained in each model, σ ∆ x , has been fitted using the expres-sion σ ∆ x / H = C σ ( ∆ x ) √ t (20)where H is the local disc thickness, the time t is measured in units c (cid:13) , 1–18 Gressel, Nelson & Turner of local orbits, and the coe ffi cients C σ ( ∆ x ) are tabulated in Tab. 2.After 200 orbits, the radial di ff usion given by (20) corresponds toa change of semimajor axis ∆ a = . ∆ a = . ∆ a = . ∆ a = . ∆ a = . ∆ a = . ff usionto be very satisfactory.Examining the longer term evolution, we use (20) to calcu-late the level of di ff usion expected over an assumed disc life timeof 5 Myr. For the fully active model, A1, we obtain ∆ a = . ∆ a = .
26 au, and for model D2 we obtain ∆ a = .
40 au. It is clear from these values that in a vertically strat-ified disc sustaining MHD turbulence without a dead zone, radialdi ff usion is predicted to be too large to be consistent with solarsystem constraints. For a disc with a significant dead zone whoseheight extends either side of the midplane by a distance ∼ H , how-ever, we see that the degree of radial mixing is strongly diminished,and the results are consistent with the distribution of asteroidal tax-onomic types (Gradie & Tedesco 1982). We now attempt to summarise our results and present a coherentpicture of how the dynamics and growth of planetesimals are af-fected by turbulent stirring in discs with and without dead zones.We frame our discussion in the context of the three planetesi-mal formation models discussed in the introduction: incremen-tal growth through particle sticking over time scales exceeding10 yr (Weidenschilling 2000; Brauer et al. 2008); concentration ofchondrules in turbulent eddies, followed by gravitational contrac-tion on the chondrule settling time – typically ∼ - 10 or-bits (Cuzzi et al. 2008); gravitational collapse of dense clumps ofmetre-sized bodies formed in turbulent discs through a combina-tion of trapping in local pressure maxima and the streaming insta-bility (Johansen et al. 2007). Although there are significant prob-lems with the incremental growth picture, in this paper we are inter-ested in exploring the e ff ects of turbulence on macroscopic bodiesof sizes R p >
100 m, and so it is useful for our discussion to assumean optimistic outcome for this model. Toward the end of this sec-tion we also discuss some of the shortcomings of our model, andissues that are raised by these for the interpretation of our results.
The results presented in Section 5 suggest that planetesimal for-mation and growth via a process involving mutual collisions be-tween smaller bodies is not possible in fully turbulent discs. Therapid excitation of large random velocities which exceed the dis-ruption / erosion threshold for planetesimals with R p
100 km willsimply lead to the destruction of bodies which form in this manner.Although the turbulence simulated in this paper is quite vigorous, because of the imposed vertical magnetic field, the scaling devel-oped in paper I for the strength of stochastic forcing as a function ofthe e ff ective viscous stress suggest that even an order of magnitudedecrease in the e ff ective viscous α value will not decrease the ran-dom velocities su ffi ciently to prevent catastrophic disruption fromoccurring.The formation of large ( R p ≃
100 km) planetesimals throughchondrules concentrating in turbulent eddies may be possible infully turbulent discs driven by the MRI. The obvious requirementfor there to be a turbulent cascade resulting in Kolmogorov-scaleeddies in which chondrules can concentrate would appear to besatisfied in such a disc. During the initial formation and settlingstage of these objects, they are likely to be of low density and hencestrongly coupled to the gas via gas drag. Relative velocities on lo-cal scales relevant for collisions are likely to be small. As they con-tract to form “sandpile” planetesimals, however, they will decouplefrom the gas and evolve dynamically like planetesimals with inter-nal densities of ρ p ≃ − . If formation / contraction times forthese bodies are ≃ yr at 5 au, then turbulent stirring will causetheir velocity dispersion to grow to the surface escape velocity of ≃
100 m s − within 3 . × yr, preventing runaway growth fromensuing. Further collisional growth between sandpile planetesimalswill therefore occur slowly. Continued driving of the velocity dis-persion by turbulence eventually allows v disp to exceed the disrup-tion value for 100 km rubble piles within approximately 1 . R p ≃
500 km). Such objects have surface escapevelocities of 500 m s − , and the time over which turbulent stirringgenerates a velocity dispersion of this magnitude is ∼ × yr. Itis unclear whether a population of planetesimals with an approxi-mately unimodal size distribution centred on R p =
500 km can un-dergo runaway growth, because of self-stirring and weak gas dragdamping of eccentricities. But if it is feasible then turbulent stirringin a fully active disc such as computed in model A1 will probablynot provide significant hindrance because of the long time scales re-quired for the velocity dispersion to grow above v disp ≃
500 m s − . Our results for model D1 show that in a model with a relativelydeep dead zone, the equilibrium velocity dispersion is determinedby collisional rather than gas drag damping for bodies with R p
10 km, for reasons discussed in section 5.3. The equilibrium ve-locity dispersion lies between the disruption / erosion thresholds forweak and strong aggregates, suggesting that incremental collisionalgrowth of planetesimals is possible, but depends on the materialstrength of the bodies. For 100 m bodies, the time required to ex- c (cid:13) , 1–18 ynamics of planetesimals in dead zones cite the velocity dispersion to the equilibrium value of 4 . − is ≃ . × yr, comfortably shorter than the formation time of km-sized planetesimals at 5 au. Similarly, the time required to excite1 km bodies to their equilibrium random velocities is ≃ . × yr,comparable to our assumed formation time for these bodies throughincremental growth. At each stage during the incremental growthof planetesimals within the dead zone, the planetesimals will havevelocity dispersions close to the equilibrium values, which thencontrol in-part the formation and growth rates. But provided thematerial strength of the planetesimals is su ffi cient (larger than thatfor weak aggregates), planetesimal growth should be possible in adead zone through collisional agglomeration, so long as the bounc-ing and metre-size barriers can be overcome (Brauer et al. 2008;Zsom et al. 2010).Runaway growth of km-sized planetesimals requires the ve-locity dispersion to be significantly lower than the escape veloc-ity from the accreting bodies, and for 10 km planetesimals v esc =
10 m s − . Results from model D1 indicate that a velocity disper-sion of this magnitude will be excited within ≃ yr, e ff ectivelyquenching the runaway growth. Continued collisional growth canstill proceed, since the velocity dispersion remains below the dis-ruption / erosion threshold (at least for strong aggregates and rubblepiles), and the equilibrium velocity dispersion for 10 km objects v disp =
51 m s − is not reached until after 2 . R p ≃
100 kmwill allow runway growth of these planetesimals to occur. ModelD1 shows that the time required to raise the velocity dispersionup to the escape velocity v esc =
100 m s − for these bodies is10 Myr. In this general picture, the influence of stochastic gravi-tational forces acting on planetesimals in a dead zone can cause thepostponement of runaway growth for 10 km planetesimals, but once100 km sized bodies form, rapid runaway growth produce numer-ous ∼ We now consider the possible e ff ects of assumptions we have madein our models that may a ff ect their outcomes and the conclusionswe have drawn from them. The disc model we adopt in this study is lower in mass by a factorof ∼ ff ects ofplanetary migration which can enable planetary embryos to growbeyond their isolation masses, thus assisting in the formation ofgiant planets within expected disc life times (Alibert et al. 2005).If relative density perturbations in the dead zone remain the samewhen the disc mass is increased, then a higher disc mass will leadto a corresponding increase in turbulent stirring, making planetes-imal destruction / erosion more likely. But, in a disc where the ion-isation sources are external (stellar X-rays, galactic cosmic rays),the column density and mass of the active zone are almost indepen-dent of the total disc mass / surface density, and depend primarilyon the penetration depth of the ionising sources. Consequently, wemay expect the relative density fluctuations near the midplane in thedead zone to be weaker in a heavier disc than in a lighter disc, po-tentially providing a less destructive source of gravitational stirring.At present it is unclear how turbulent stirring in a dead zone scaleswith the disc mass, and studying this issue is di ffi cult from the com-putational point of view. The shearing box models we present hereare extended in the x and y directions compared to those used inmost studies, and include ± . We have assumed a particular configuration for the initial magneticfields in our simulations that includes a small but non-negligiblenet-flux vertical component. It is well known that the existenceand strength of net-flux magnetic fields changes the saturationlevel of developed MHD turbulence in simulations of the MRI(Hawley et al. 1995). A future study of planetesimal dynamics em-bedded in turbulent discs with dead zones should examine the roleplayed by the field topology and strength.
In our estimates of collisional damping we assume that all plan-etesimals are the same size, and that all solids are incorporated intoplanetesimals. The latter assumption is broadly consistent with ourdisc model, in which we assume that 90% of the dust has beendepleted due to the growth of grains into larger bodies. But theassumption of a single planetesimal size is not realistic. A more ac-curate calculation of the role of collisional damping would requirea self-consistent model of accretion and fragmentation such that the c (cid:13) , 1–18 Gressel, Nelson & Turner size distribution is modelled self-consistently. Unfortunately sucha model goes beyond the scope of this paper.
When estimating the scale height of the planetesimal disc dur-ing the calculation of collisional damping rates, we assumed that H p = v disp / Ω k (i.e. the velocity dispersion is isotropic). In a turbu-lent disc, however, where a significant component of the stochasticgravitational field is generated by spiral density waves that have lit-tle vertical structure, we should expect that the eccentricity growthexceeds that of the inclination. Our simulation results, however,show that for models D1 and D2 the inclination growth actuallyexceeds that of the eccentricity. This surprising result arises be-cause the disc midplane (centre of mass in the z -direction), oscil-lates vertically with an amplitude that reaches ∼ . H , causing anoscillation of similar magnitude in the planetesimals. Future simu-lations of planetesimals within dead zones will be computed usingglobal models to examine how robust this phenomenon is, althoughwe note that a similar vertical oscillation has been reported in theshearing box simulations presented by Oishi et al. (2007). Wave propagation in discs is strongly influenced by the equation ofstate (Lin et al. 1990; Korycansky & Pringle 1995). This may a ff ectthe ability of waves that are excited in the active regions to propa-gate into the dead zone, and thus a ff ect the turbulent stirring there.The simulations presented in this paper all adopted an isothermalequation of state, whereas protoplanetary discs are expected to beoptically thick at 5 au.For parameters typical of T Tauri stars, the primary heatingmechanism near the midplane is the absorption of the radiationemitted by the disc atmosphere, which in turn heats mainly byabsorbing light from the central star (Chiang & Goldreich 1997).The accretion heating also occurs in the atmosphere, leading toa thoroughly isothermal interior when a dead zone is present(Hirose & Turner 2011). The role of the equation of state in de-termining the level of turbulent stirring at the midplane will be ex-plored in a future publication. In this paper we have presented the results of simulations of plan-etesimals embedded in local models of protoplanetary discs. Themain aim of this work is to understand the dynamics of planetesi-mals in turbulent discs, and to examine how dead zones of varyingvertical thicknesses change the stochastic forcing of the radial ve-locity dispersion, and the di ff usion of planetesimal semimajor axes.We consider the influence of gas drag and collisional damping incounteracting the growth of random velocities due to turbulent stir-ring, and estimate the equilibrium velocity dispersions which re-sult. The main conclusions of this study are as follows.The presence and depth of a dead zone have negligible e ff ecton the measured coherence time of the particle stirring, but onlya ff ect its amplitude as measured through the width of the torquedistribution.In fully turbulent discs without dead zones, we find that thevelocity dispersion of embedded planetesimals grows rapidly, andquickly exceeds the threshold for disruption for planetesimals in the size range 100 m R p
10 km. We conclude that planetesimal for-mation via collisional accretion of smaller bodies cannot occur inglobally turbulent discs. The direct formation of large 100 km-sizedplanetesimals via turbulent concentration and gravitational insta-bility (Johansen et al. 2007; Cuzzi et al. 2008) may occur, but theexcitation of their random velocities may prevent runaway growth,making it di ffi cult to form large planetary embryos and cores withinreasonable disc life times.Radial di ff usion of planetesimals occurs over distances of ∼ . . ± H above and below the midplane, the stochastic ex-citation of planetesimal eccentricities is reduced by a factor ∼ R p =
10 km bodies within10 yr, quenching runaway growth. Slower growth of planetesimalsthrough collisions remains possible, and runaway growth can ensueonce bodies with R p ≃
100 km form, due to the slow growth of ran-dom velocities within a dead zone. If R p ≃ • Comparing the results obtained for dead zones with di ff erentvertical heights ( ∼ . H and 2 H ), we find that the larger dead zonemodel clearly favours the collisional formation of planetesimals be-cause of weaker stochastic forcing of random motions. • The radial di ff usion of planetesimals is much reduced in discswith dead zones relative to fully active discs. Di ff usion over a dis-tance equal to 0 .
26 au is predicted by our deeper dead zone model,consistent with solar system constraints. • The low di ff usion levels for planetesimal semimajor axes indead zones clearly means that stochastic torques in such an envi-ronment cannot prevent the large scale inward drift of low massplanets via type I migration. c (cid:13) , 1–18 ynamics of planetesimals in dead zones • We suggest that dead zones provide safe havens for km-sizedplanetesimals against destructive collisions, and are probably a re-quired ingredient for the formation of planetary systems.This paper is part of series in which we are examining thedynamics of planetesimals and planets in increasingly realistic discmodels, and shows that di ff erent qualitative outcomes for planetformation can be expected in discs with dead zones versus thosewithout. In future publications we will examine a number of issues,including how changing the disc surface density and the mass ratiobetween the active and dead zones modifies the turbulent stirring,and how turbulence varies as a function of heliocentric distance inglobal MHD simulations of protoplanetary discs. ACKNOWLEDGEMENTS
This work used the nirvana - iii code developed by Udo Ziegler at theAstrophysical Institute Potsdam. All computations were performedon the QMUL HPC facility, purchased under the SRIF initiative.R.P.N. and O.G. acknowledge the hospitality of the Isaac New-ton Institute for Mathematical Sciences, where part of the workpresented in this paper was completed during the ‘Dynamics ofDiscs and Planets’ research programme. N.J.T. was supported bythe Jet Propulsion Laboratory, California Institute of Technology,the NASA Origins and Outer Planets programs, and the Alexan-der von Humboldt Foundation. We thank the referee, Stuart Wei-denschilling, for useful comments that led to improvements to thispaper. REFERENCES
Alibert Y., Mordasini C., Benz W., Winisdoer ff er C., 2005, A&A,434, 343Armitage P. J., 1998, ApJ, 501, L189 + Asphaug E., Benz W., 1996, Icarus, 121, 225Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Balsara D. S., Meyer C., 2010, (astro-ph:1003.0018)Baruteau C., Lin D. N. C., 2010, ApJ, 709, 759Benz W., Asphaug E., 1999, Icarus, 142, 5Blaes O. M., Balbus S. A., 1994, ApJ, 421, 163Brandenburg A., 2008, AN, 329, 725Brandenburg A., Candelaresi S., Chatterjee P., 2009, MNRAS,398, 1414Brauer F., Dullemond C. P., Henning T., 2008, A&A, 480, 859Bridges F. G., Hatzes A., Lin D. N. C., 1984, Nature, 309, 333Chiang E. I., Goldreich P., 1997, ApJ, 490, 368Courvoisier A., Hughes D. W., Proctor M. R. E., 2010, AN, 331,667Cuzzi J. N., Hogan R. C., Shari ff K., 2008, ApJ, 687, 1432Davis S. W., Stone J. M., Pessah M. E., 2010, ApJ, 713, 52Flaig M., Kley W., Kissmann R., 2010, MNRAS, 409, 1297Fleming T., Stone J. M., 2003, ApJ, 585, 908Fleming T. P., Stone J. M., Hawley J. F., 2000, ApJ, 530, 464Flock M., Dzyurkevich N., Klahr H., Mignone A., 2010, A&A,516, A26 + Fromang S., Terquem C., Balbus S. A., 2002, MNRAS, 329, 18Gammie C. F., 1996, ApJ, 457, 355Gardiner T. A., Stone J. M., 2005, JCoPh, 205, 509Garmire G., Feigelson E. D., Broos P., Hillenbrand L. A., PravdoS. H., Townsley L., Tsuboi Y., 2000, AJ, 120, 1426Glassgold A. E., Najita J., Igea J., 1997, ApJ, 480, 344 Goldreich P., Ward W. R., 1973, ApJ, 183, 1051Gradie J., Tedesco E., 1982, Science, 216, 1405Gressel O., 2010, MNRAS, 405, 41Gressel O., Ziegler U., 2007, CoPhC, 176, 652Haisch Jr. K. E., Lada E. A., Lada C. J., 2001, ApJ, 553, L153Hanasz M., Otmianowska-Mazur K., Kowal G., Lesch H., 2009,A&A, 498, 335Hawley J. F., 2001, ApJ, 554, 534Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742Hayashi C., 1981, Progress of Th. Phys. Suppl., 70, 35Heinemann T., Papaloizou J. C. B., 2009a, MNRAS, 397, 52Heinemann T., Papaloizou J. C. B., 2009b, MNRAS, 397, 64Hirose S., Turner N. J., 2011, ApJ, 732, L30Ida S., Guillot T., Morbidelli A., 2008, ApJ, 686, 1292Ida S., Makino J., 1993, Icarus, 106, 210Igea J., Glassgold A. E., 1999, ApJ, 518, 848Ilgner M., Nelson R. P., 2006a, A&A, 445, 205Ilgner M., Nelson R. P., 2006b, A&A, 445, 223Ilgner M., Nelson R. P., 2008, A&A, 483, 815Johansen A., Klahr H., Henning T., 2011, A&A, 529, A62 + Johansen A., Levin Y., 2008, A&A, 490, 501Johansen A., Oishi J. S., Low M., Klahr H., Henning T., YoudinA., 2007, Nature, 448, 1022Johansen A., Youdin A., Mac Low M., 2009, ApJ, 704, L75Johnson E. T., Goodman J., Menou K., 2006, ApJ, 647, 1413Kenyon S. J., Bromley B. C., 2009, ApJ, 690, L140Kokubo E., Ida S., 1998, Icarus, 131, 171Kortenkamp S. J., Wetherill G. W., Inaba S., 2001, Science, 293,1127Korycansky D. G., Pringle J. E., 1995, MNRAS, 272, 618Latter H. N., Fromang S., Gressel O., 2010, MNRAS, 406, 848Latter H. N., Lesa ff re P., Balbus S. A., 2009, MNRAS, 394, 715Lin D. N. C., Papaloizou J. C. B., Savonije G. J., 1990, ApJ, 364,326Miller K. A., Stone J. M., 2000, ApJ, 534, 398Miyoshi T., Kusano K., 2005, JCoPh, 208, 315Moth´e-Diniz T., Carvano J. M. ´A., Lazzaro D., 2003, Icarus, 162,10Nelson R. P., 2005, A&A, 443, 1067Nelson R. P., Gressel O., 2010, MNRAS, 409, 639Nelson R. P., Papaloizou J. C. B., 2004, MNRAS, 350, 849O’Brien D. P., Morbidelli A., Levison H. F., 2006, Icarus, 184, 39Oishi J. S., Mac Low M., Menou K., 2007, ApJ, 670, 805Papaloizou J. C. B., Nelson R. P., 2003, MNRAS, 339, 983Pessah M. E., Goodman J., 2009, ApJ, 698, L72Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J.,Podolak M., Greenzweig Y., 1996, Icarus, 124, 62Regev O., Umurhan O. M., 2008, A&A, 481, 21Rheinhardt M., Brandenburg A., 2010, A&A, 520, A28 + Richardson J. E., Melosh H. J., Lisse C. M., Carcich B., 2007,Icarus, 190, 357R¨udiger G., Pipin V. V., 2000, A&A, 362, 756Sano T., Miyama S. M., Umebayashi T., Nakano T., 2000, ApJ,543, 486Shi J., Krolik J. H., Hirose S., 2010, ApJ, 708, 1716Sicilia-Aguilar A., Hartmann L. W., Brice˜no C., Muzerolle J., Cal-vet N., 2004, AJ, 128, 805Stewart S. T., Leinhardt Z. M., 2009, ApJ, 691, L133Stone J. M., Gardiner T. A., 2010, ApJS, 189, 142Stone J. M., Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ,463, 656Turner N. J., Drake J. F., 2009, ApJ, 703, 2152 c (cid:13) , 1–18 Gressel, Nelson & Turner
Turner N. J., Sano T., 2008, ApJ, 679, L131Turner N. J., Sano T., Dziourkevitch N., 2007, ApJ, 659, 729Umebayashi T., 1983, Progress of Theoretical Physics, 69, 480Umebayashi T., Nakano T., 1981, PASJ, 33, 617Umebayashi T., Nakano T., 2009, ApJ, 690, 69Wardle M., 1999, MNRAS, 307, 849Weidenschilling S. J., 1977, MNRAS, 180, 57Weidenschilling S. J., 2000, Space Sci. Rev., 92, 295Weidenschilling S. J., Cuzzi J. N., 1993, in E. H. Levy & J. I. Lu-nine eds., Protostars and Planets III, Univ. Arizona Press, Tuc-son, p. 1031Wetherill G. W., Stewart G. R., 1993, Icarus, 106, 190Yang C., Mac Low M., Menou K., 2009, ApJ, 707, 1233Youdin A. N., Goodman J., 2005, ApJ, 620, 459Ziegler U., 2004, JCoPh, 196, 393Ziegler U., 2008, CoPhC, 179, 227Zsom A., Ormel C. W., G¨uttler C., Blum J., Dullemond C. P.,2010, A&A, 513, A57 + This paper has been typeset from a TEX / L A TEX file prepared by theauthor. c (cid:13)000