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 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 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&

On the ionisation fraction in protoplanetary disks III. The effect of X-ray flares on gas-phase chemistry

Context. Recent observations of the X-ray emission from T Tauri stars in the Orion nebula have shown that they undergo frequent outbursts in their X-ray luminosity. These X-ray flares are characterised by increases in luminosity by two orders of magnitude, a typical duration of less than one day, and a significant hardening of the X-ray spectrum. Aims. It is unknown what effect these X-ray flares will have on the ionisation fraction and dead-zone structure in protoplanetary disks. We present the results of calculations designed to address this question. Methods. We have performed calculations of the ionisation fraction in a standard $\alpha$-disk model using two different chemical reaction networks. We include in our models ionisation due to X-rays from the central star, and calculate the time-dependent ionisation fraction and dead--zone structure for the inner 10 AU of a protoplanetary disk model. Results. We find that the disk response to X-ray flares depends on whether the plasma temperature increases during flares and/or whether heavy metals (such as magnesium) are present in the gas phase. Under favourable conditions the outer disk dead--zone can disappear altogether,and the dead-zone located between 0.5 < R < 2 AU can disappear and reappear in phase with the X-ray luminosity. Conclusions. X-ray flares can have a significant effect on the dead-zone structure in protoplanetary disks. Caution is required in interpreting this result as the duration of X-ray bursts is considerably shorter than the growth time of MHD turbulence due to the magnetorotational instability.Comment: 12 pages, 8 figures, accepted by A &

On the growth and orbital evolution of giant planets in layered protoplanetary disks

We present the results of hydrodynamic simulations of the growth and orbital evolution of giant planets embedded in a protoplanetary disk with a dead-zone. The aim is to examine to what extent the presence of a dead-zone affects the rates of mass accretion and migration for giant planets. We performed 3D numerical simulations using a grid-based hydrodynamics code. In these simulations of non-magnetised disks, the dead-zone is treated as a region where the vertical profile of the viscosity depends on the distance from the equatorial plane. We consider dead-zones with vertical sizes, H_dz, ranging from 0 to H_dz=2.3H, where H is the disk scale-height. For all models, the vertically integrated viscous stress, and the related mass flux through the disk, have the same value, such that the simulations test the dependence of planetary mass accretion and migration on the vertical distribution of the viscous stress. For each model, an embedded 30 earth-masses planet on a fixed circular orbit is allowed to accrete gas from the disk. Once the planet mass becomes equal to that of Saturn or Jupiter, we allow the planet orbit to evolve due to gravitational interaction with the disk. We find that the time scale over which a protoplanet grows to become a giant planet is essentially independent of the dead-zone size, and depends only on the total rate at which the disk viscously supplies material to the planet. For Saturn-mass planets, the migration rate depends only weakly on the size of the dead-zone for H_dz< 1.5H, but becomes slower when H_dz=2.3H. This effect is due to the desaturation of corotation torques which originate from residual material in the partial-gap region. For Jupiter-mass planets, there is a clear tendency for the migration to proceed more slowly as the size of the dead-zone increases.Comment: Accepted for publication in A&A. 10 pages, 12 figure