Global MHD simulations of stratified and turbulent protoplanetary discs. I. Model properties

We present the results of global 3-D MHD simulations of stratified and turbulent protoplanetary disc models. The aim of this work is to develop thin disc models capable of sustaining turbulence for long run times, which can be used for on-going studies of planet formation in turbulent discs. The results are obtained using two codes written in spherical coordinates: GLOBAL and NIRVANA. Both are time--explicit and use finite differences along with the Constrained Transport algorithm to evolve the equations of MHD. In the presence of a weak toroidal magnetic field, a thin protoplanetary disc in hydrostatic equilibrium is destabilised by the magnetorotational instability (MRI). When the resolution is large enough (25 vertical grid cells per scale height), the entire disc settles into a turbulent quasi steady-state after about 300 orbits. Angular momentum is transported outward such that the standard alpha parameter is roughly 4-6*10^{-3}. We find that the initial toroidal flux is expelled from the disc midplane and that the disc behaves essentially as a quasi-zero net flux disc for the remainder of the simulation. As in previous studies, the disc develops a dual structure composed of an MRI--driven turbulent core around its midplane, and a magnetised corona stable to the MRI near its surface. By varying disc parameters and boundary conditions, we show that these basic properties of the models are robust. The high resolution disc models we present in this paper achieve a quasi--steady state and sustain turbulence for hundreds of orbits. As such, they are ideally suited to the study of outstanding problems in planet formation such as disc--planet interactions and dust dynamics.Comment: 19 pages, 29 figures, accepted in Astronomy & Astrophysic

The interaction of a giant planet with a disc with MHD turbulence I: The initial turbulent disc models

This is the first of a series of papers aimed at developing and interpreting simulations of protoplanets interacting with turbulent accretion discs. Here we study the disc models prior to the introduction of a protoplanet.We study models in which a Keplerian domain is unstable to the magnetorotational instability (MRI). Various models with B-fields having zero net flux are considered.We relate the properties of the models to classical viscous disc theory.All models attain a turbulent state with volume averaged stress parameter alpha ~ 0.005. At any particular time the vertically and azimuthally averaged value exhibited large fluctuations in radius. Time averaging over periods exceeding 3 orbital periods at the outer boundary of the disc resulted in a smoother quantity with radial variations within a factor of two or so. The vertically and azimuthally averaged radial velocity showed much larger spatial and temporal fluctuations, requiring additional time averaging for 7-8 orbital periods at the outer boundary to limit them. Comparison with the value derived from the averaged stress using viscous disc theory yielded schematic agreement for feasible averaging times but with some indication that the effects of residual fluctuations remained. The behaviour described above must be borne in mind when considering laminar disc simulations with anomalous Navier--Stokes viscosity. This is because the operation of a viscosity as in classical viscous disc theory with anomalous viscosity coefficient cannot apply to a turbulent disc undergoing rapid changes due to external perturbation. The classical theory can only be used to describe the time averaged behaviour of the parts of the disc that are in a statistically steady condition for long enough for appropriate averaging to be carried out.Comment: 10 pages, 23 figures, accepted for publication in MNRAS. A gzipped postscript version including high resolution figures is available at http://www.maths.qmul.ac.uk/~rp

On the Ionisation Fraction in Protoplanetary Disks I: Comparing Different Reaction Networks

We calculate the ionisation fraction in protostellar disk models using a number of different chemical reaction networks, including gas-phase and gas-grain reaction schemes. The disk models we consider are conventional alpha-disks, which include viscous heating and radiative cooling. The primary source of ionisation is assumed to be X-ray irradiation from the central star. We consider a number of gas-phase chemical networks. In general we find that the simple models predict higher fractional ionisation levels and more extensive active zones than the more complex models. When heavy metal atoms are included the simple models predict that the disk is magnetically active throughout. The complex models predict that extensive regions of the disk remain magnetically uncoupled even with a fractional abundance of magnesium of 10(-8). The addition of submicron sized grains with a concentration of 10(-12) causes the size of the dead zone to increase dramatically for all kinetic models considered. We find that the simple and complex gas-grain reaction schemes agree on the size and structure of the resulting dead zone. We examine the effects of depleting the concentration of small grains as a crude means of modeling the growth of grains during planet formation. We find that a depletion factor of 10(-4) causes the gas-grain chemistry to converge to the gas-phase chemistry when heavy metals are absent. 10(-8) is required when magnesium is included. This suggests that efficient grain growth and settling will be required in protoplanetary disks, before a substantial fraction of the disk mass in the planet forming zone between 1 - 10 AU becomes magnetically active and turbulent.Comment: 21 pages, 23 figures, accepted for publication in A & A Includes correction to our implementation of the Umebayashi-Nakano reaction networ

On the Orbital Evolution of Low Mass Protoplanets in Turbulent, Magnetised Disks

(Abridged).We present the results of MHD simulations of low mass protoplanets interacting with turbulent disks. We calculate the orbital evolution of `planetesimals' and protoplanets with masses in the range 0 < m_p < 30 M_Earth. Planetesimals and protoplanets undergo stochastic migration due to interaction with turbulent density fluctuations. Over run times of ~ 150 planet orbits, stochastic migration dominates over type I migration for many models. Fourier analysis of the torques experienced by planets indicates that the torque fluctuations contain components with significant power whose time scales of variation are similar to the simulation run times. These low frequency fluctuations partly explain the dominance of stochastic torques, and may provide a powerful means of counteracting the type I migration of some planets in turbulent disks. Turbulence is a source of eccentricity driving. Planetesimals attained eccentricities in the range 0.02 < e < 0.14, m_p=1 M_Earth planets attained eccentricities 0.02 < e < 0.08, and m_p=10 M_Earth protoplanets reached 0.02 < e < 0.03. This is in basic agreement with a model in which turbulence drives e-growth, and interaction with disk material at coorbital Lindblad resonances causes e-damping. These results are significant for planet formation. Stochastic migration may prevent some planet cores migrating into their star via type I before becoming gas giants. The growth of planetary cores may be enhanced by preventing isolation. Eccentricity excitation by turbulence, however, may reduce growth rates of planetary cores during the runaway and oligarchic growth stages, and cause collisions between planetesimals to become destructive.Comment: 21 pages, 16 figures. Accepted for publication in Astronomy & Astrophysics. A version with full resolution, colour figures is available from: http://www.maths.qmul.ac.uk/~rpn/preprint

The interaction of planets with a disc with MHD turbulence III: Flow morphology and conditions for gap formation in local and global simulations

We present the results of both global cylindrical disc simulations and local shearing box simulations of protoplanets interacting with a disc undergoing MHD turbulence with zero net flux magnetic fields. We investigate the nature of the disc response and conditions for gap formation. This issue is an important one for determining the type and nature of the migration of the protoplanet, with the presence of a deep gap being believed to enable slower migration. For both types of simulation we find a common pattern of behaviour for which the main parameter determining the nature of the response is $M_p R^3/(M_* H^3)$, with $M_p$, $M_*$, $R$, and $H$ being the protoplanet mass, the central mass, the orbital radius and the disc semi-thickness respectively. We find that as this parameter is increased towards 0.1, the presence of the protoplanet is first indicated by the appearance of the well known trailing wake which, although it may appear erratic on account of the turbulence, appears to be well defined. Once the above parameter exceeds a number around unity a gap starts to develop inside which the magnetic energy density tends to be concentrated in the high density wakes. This gap formation condition can be understood from simple dimensional considerations of the conditions for nonlinearity, and the balance of angular momentum transport due to Maxwell and Reynolds' stresses with that due to tidal torques. An important result is that the basic flow morphology in the vicinity of the protoplanet is very similar in both the local and global simulations. This indicates that local shearing box simulations, which are computationally less demanding, capture much of the physics of disc-planet interaction. Thus they may provide a useful tool for studying the local interaction between forming protoplanets and turbulent, protostellar discs.Comment: 20 pages, 28 figures (some colour), accepted for publication in M.N.R.A.S. with minor modification. A pdf version containing high resolution colour figures is available from http://www.maths.qmul.ac.uk/~rpn/projects/mhd along with additional images and movies. A companion paper accepted without change by M.N.R.A.S. is also availabl

The interaction of a giant planet with a disc with MHD turbulence II: The interaction of the planet with the disc

We present a global MHD simulation of a turbulent accretion disc interacting with a protoplanet of 5 Jupiter masses. The disc model had H/r=0.1,and a value of the Shakura & Sunyaev alpha ~ 0.005. The protoplanet opened a gap in the disc, with the interaction leading to inward migration on the expected time scale. Spiral waves were launched by the protoplanet and although they were diffused and dissipated through interaction with the turbulence, they produced an outward angular momentum flow which compensated for a reduced flux associated with the turbulence, so maintaining the gap. When compared with laminar disc models with the same estimated alpha, the gap was found to be deeper and wider indicating that the turbulent disc behaved as if it possessed a smaller alpha. This may arise for two reasons. First, the turbulence does not provide a source of constantly acting friction in the near vicinity of the planet that leads to steady mass flow into the gap region. Instead the turbulence is characterised by large fluctuations in the radial velocity, and time averaging over significant time scales is required to recover the underlying mass flow through the disc. Near the planet the disc material experiences high amplitude perturbations on time scales that are short relative to the time scale required for averaging. The disc response is thus likely to be altered relative to a Navier--Stokes model. Second, the simulation indicates that an ordered magnetic connection between the inner and outer disc can occur enabling angular momentum to flow out across the gap, helping to maintain it independently of the protoplanet's tide. This type of effect may assist gap formation for smaller mass protoplanets which otherwise would not be able to maintain them.Comment: 14 pages, 17 figures, accepted for publication in MNRAS. A gzipped postscript version including high resolution figures is available at http://www.maths.qmw.ac.uk/~rp

On the Ionisation Fraction in Protoplanetary Disks II: The Effect of Turbulent Mixing on Gas--phase Chemistry

We calculate the ionisation fraction in protostellar disk models using two different gas-phase chemical networks, and examine the effect of turbulent mixing by modelling the diffusion of chemical species vertically through the disk. The aim is to determine in which regions of the disk gas can couple to a magnetic field and sustain MHD turbulence. We find that the effect of diffusion depends crucially on the elemental abundance of heavy metals (magnesium) included in the chemical model. In the absence of heavy metals, diffusion has essentially no effect on the ionisation structure of the disks, as the recombination time scale is much shorter than the turbulent diffusion time scale. When metals are included with an elemental abundance above a threshold value, the diffusion can dramatically reduce the size of the magnetically decoupled region, or even remove it altogther. For a complex chemistry the elemental abundance of magnesium required to remove the dead zone is 10(-10) - 10(-8). We also find that diffusion can modify the reaction pathways, giving rise to dominant species when diffusion is switched on that are minor species when diffusion is absent. This suggests that there may be chemical signatures of diffusive mixing that could be used to indirectly detect turbulent activity in protoplanetary disks. We find examples of models in which the dead zone in the outer disk region is rendered deeper when diffusion is switched on. Overall these results suggest that global MHD turbulence in protoplanetary disks may be self-sustaining under favourable circumstances, as turbulent mixing can help maintain the ionisation fraction above that necessary to ensure good coupling between the gas and magnetic field.Comment: 11 pages, 7 figures; accepted for publication in A &

The interaction of planets with a disc with MHD turbulence IV: Migration rates of embedded protoplanets

(Abridged) We present global disc and local shearing box simulations of planets interacting with a MHD turbulent disc. We examine the torque exerted by the disc on the embedded planets as a function of planet mass, and thus make a first study of orbital migration of planets due to interaction with turbulent discs. Global simulations were performed for a disc with H/R=0.07 and planet masses M_p=3,10,30 Earth masses, and 3 Jupiter masses. Shearing box runs were performed for values of (M_p/M_*)/(H/R)^3=0.1,0.3,1.0 and 2.0, M_* being the central mass. These allow embedded and gap forming planets to be examined. In all cases the instantaneous torque exerted on a planet showed strong fluctuations. In the embedded cases it oscillated between negative and positive values, and migration occurs as a random walk, unlike the usual type I migration. Running time averages for embedded planets over 20-25 orbital periods show that large fluctuations occur on longer time scales, preventing convergence of the average torque to well defined values, or even to a well defined sign. Fluctuations become relatively smaller for larger masses, giving better convergence, due to the planet's perturbation of the disc becoming larger than the turbulence in its vicinity. Eventually gap formation occurs, with a transition to type II migration. The existence of significant fluctuations occurring in turbulent discs on long time scales is important for lower mass embedded protoplanets. If significant fluctuations occur on the longest disc evolutionary time scales, convergence of torque running averages for practical purposes will not occur, and the migration behaviour of low mass protoplanets considered as an ensemble would be very different from predictions of type I theory for laminar discs.Comment: 19 pages, 24 figures (some colour), submitted to M.N.R.A.S. A gzipped postscript version containing high resolution colour figures is available from http://www.maths.qmul.ac.uk/~rp

On the dynamics of planetesimals embedded in turbulent protoplanetary discs

(abridged) Angular momentum transport and accretion in protoplanetary discs are generally believed to be driven by MHD turbulence via the magneto-rotational instability (MRI). The dynamics of solid bodies embedded in such discs (dust grains, boulders, planetesimals and planets) may be strongly affected by the turbulence, such that the formation pathways for planetary systems are determined in part by the strength and spatial distribution of the turbulent flow. We examine the dynamics of planetesimals, with radii between 1m \^a 10 km, embedded in turbulent protoplanetary discs, using three dimensional MHD simulations. The planetesimals experience gas drag and stochastic gravitational forces due to the turbulent disc. We use, and compare the results from, local shearing box simulations and global models in this study. The main aims of this work are to examine: the growth, and possible saturation, of the velocity dispersion of embedded planetesimals as a function of their size and disc parameters; the rate of radial migration and diffusion of planetesimals; the conditions under which the results from shearing box and global simulations agree. We find good agreement between local and global simulations when shearing boxes of dimension 4H x 16H x 2H are used (H being the local scale height). The magnitude of the density fluctuations obtained is sensitive to the box size, due to the excitation and propagation of spiral density waves. This affects the stochastic forcing experienced by planetesimals. [...] Our models show that fully developed MHD turbulence in protoplanetary discs would have a destructive effect on embedded planetesimals. Relatively low levels of turbulence are required for traditional models of planetesimal accretion to operate, this being consistent with the existence of a dead zone in protoplanetary discs.Comment: 23 pages, 28 figures, 3 tables, accepted for publication in MNRA

On the migration of protoplanets embedded in circumbinary disks

We present the results of hydrodynamical simulations of low mass protoplanets embedded in circumbinary accretion disks. The aim is to examine the migration and long term orbital evolution of the protoplanets, in order to establish the stability properties of planets that form in circumbinary disks. Simulations were performed using a grid--based hydrodynamics code. First we present a set of calculations that study how a binary interacts with a circumbinary disk. We evolve the system for 10^5 binary orbits, which is the time needed for the system to reach a quasi-equilibrium state. From this time onward the apsidal lines of the disk and the binary are aligned, and the binary eccentricity remains essentially unchanged with a value of e_b ~ 0.08. Once this stationary state is obtained, we embed a low mass protoplanet in the disk and let it evolve under the action of the binary and disk forces. We consider protoplanets with masses of 5, 10 and 20 Earth masses. In each case, we find that inward migration of the protoplanet is stopped at the edge of the tidally truncated cavity formed by the binary. This effect is due to positve corotation torques, which can counterbalance the net negative Lindblad torques in disk regions where the surface density profile has a sufficiently large positive gradient. Halting of migration occurs in a region of long-term stability, suggesting that low mass circumbinary planets may be common, and that gas giant circumbinary planets should be able to form in circumbinary disks.Comment: 10 pages, 10 figures, accepted for publication in A&