EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 497 / No 2 (April II 2009) A&A, 497 2 (2009) 595-609 Full HTML
Free Access
Issue
A&A
Volume 497, Number 2, April II 2009
Page(s) 595 - 609
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200811330
Published online 18 February 2009

On the evolution of mean motion resonances through stochastic forcing: fast and slow libration modes and the origin of HD 128311

H. Rein - J. C. B. Papaloizou

Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK

Received 11 November 2008 / Accepted 10 February 2009

Abstract
Aims. We clarify the response of extrasolar planetary systems in a 2:1 mean motion commensurability with masses ranging from the super Jovian range to the terrestrial range to stochastic forcing that could result from protoplanetary disk turbulence. The behaviour of the different libration modes for a wide range of system parameters and stochastic forcing magnitudes is investigated. The growth of libration amplitudes is parameterized as a function of the relevant physical parameters. The results are applied to provide an explanation of the configuration of the HD 128311 system.
Methods. We first develop an analytic model from first principles without making the assumption that both eccentricities are small. We also perform numerical N-body simulations with additional stochastic forcing terms to represent the effects of putative disk turbulence.
Results. We isolate two distinct libration modes for the resonant angles. These react to stochastic forcing in a different way and become coupled when the libration amplitudes are large. Systems are quickly destabilized by large magnitudes of stochastic forcing but some stability is imparted should systems undergo a net orbital migration. The slow mode, which mostly corresponds to motion of the angle between the apsidal lines of the two planets, is converted to circulation more readily than the fast mode which is associated with oscillations of the semi-major axes. This mode is also vulnerable to the attainment of small eccentricities which causes oscillations between periods of libration and circulation.
Conclusions. Stochastic forcing due to disk turbulence may have played a role in shaping the configurations of observed systems in mean motion resonance. It naturally provides a mechanism for accounting for the HD 128311 system for which the fast mode librates and the slow mode is apparently near the borderline between libration and circulation.

Key words: turbulence - celestial mechanics - planetary systems: formation - planetary systems: protoplanetary disks - methods: analytical - methods: N-body simulations

1 Introduction

Of the recently discovered 335 extrasolar planets, at least 75 are in multiple planet systems (Schneider 2009). About 10% of these are in or very close to a resonant configuration where two planets show a mean motion commensurability, with four systems in or near a 2:1 resonance (Reipurth et al. 2007; Udry et al. 2007).

Resonant configurations can be established by dissipative forces acting on the planets which lead to convergent migration (see for example Snellgrove et al. 2001; Lee & Peale 2001). An anomalous effective kinematic viscosity $\nu$ of order  ${\sim} 10^{-5}$, but with considerable uncertainty, has been inferred from observations of accretion rates onto the central protostars. The magneto rotational instability (MRI) is thought to be responsible for this anomalous value of $\nu$ which is often characterized using the Shakura & Sunyaev (1973) $\alpha =\nu/(c_{\rm s} \Omega)$ parameter, with $c_{\rm s}$ being the local sound speed and $2\pi/\Omega$ being the orbital period. For the most likely situation of small or zero net magnetic flux, recent MHD simulations have indicated $\alpha\approx 10^{-2} {-} 10^{-3}$. Even assuming adequate ionization, this value is uncertain when realistic values of the actual transport coefficients are employed because of numerical resolution issues (see Fromang & Papaloizou 2007; Fromang et al. 2007). Which parts of protoplanetary disks are adequately ionized, or consitute a dead zone, is also an issue (e.g. Gammie 1996). Because the level of MRI turbulence and associated density fluctuations are very uncertain, we consider stochastic force amplitudes ranging over several orders of magnitude.

The influence of stochastic forces resulting from the gravitational field produced by the density fluctuations associated with MRI turbulence on migrating planets was explored first by Nelson & Papaloizou (2004). They considered MRI simulations directly with the result that the simulation ran only for a relatively small number of orbits. To consider much longer evolution times, as done here, a simpler parameterized model of the forcing is needed.

Recent studies by Adams et al. (2008) and Lecoanet et al. (2008) have been made in order to estimate the lifetimes of the mean motion resonances in the (Obj)GJ 876(Obj) system when it is perturbed by a sequence of stochastic kicks. They found that resonances were disrupted within expected disk lifetimes for sufficiently large forcing. They assume a fixed orbit for the outer planet and gave an interesting discussion of the interactions of the two planets based on a pendulum equation with an additional ad hoc stochastic forcing term. The aim of the work presented in this paper is to provide a more complete study of the full set of celestial mechanics equations and to describe the interplay between apsidal corotation and mean motion resonance as well as to consider a wide range of planet masses and stochastic force amplitudes. In addition, we derive a prescription for incorporating continuous stochastic forcing terms that allows for a general autocorrelation function with associated correlation time. Their magnitudes are related to properties of the protoplanetary disk and physical scaling laws are found. This is of particular importance as the factors that control the strength of the stochastic forcing are not well constrained. We also use our formalism to estimate the stochastic diffusion rates of the orbital elements from first principles.

We further remark that Adams et al. (2008) gave a discussion of the diffusion of resonant angles that do not satisfy the d'Alambert condition (see e.g. Hamilton 1994) which is equivalent to the requirement of rotational invariance. Thus their appearrance cannot be straightforwardly connected to the basic equations or angles discussed in this paper unless one assumes that the variation of the longtitude of pericentre of one or both of the interacting planets can be neglected. However, the system lifetimes derived from their numerical work is broadly consistent with ours for the parameter regime they considered.

We find that stochastic forcing readily produces systems in mean motion resonance with broken apsidal corotation. An additional aim of this paper is to use this feature to construct scenarios involving convergent migration and stochastic forcing to account for the (Obj)HD 128311(Obj) system.

The plan of the paper is as follows. In Sect. 2 we present the basic equations governing two interacting planets subject to external stochastic forces. We specialize to the case where the planets are either in or near a mean motion commensurability and so retain only the two angle variables that vary on a time scale much longer than the orbital one. Considering the case where these angles undergo small amplitude librations, we identify fast and slow libration modes, the former being associated with variations in the semi-major axes and the latter with the angle between the apsidal lines of the two planets. We go on to consider the effects of stochastic forcing arising from some external process such as disk turbulence in Sect. 2.3 deriving the diffusion rates for the orbital elements of a single planet and the growth rate of the libration amplitudes in the two planet case.

In Sect. 3 we discuss the origin and numerical implementation of the stochastic forcing and their operation in the single planet case. In Sect. 4 we go on to consider the stochastic forcing in the two planet case. We consider the conversion of the fast and slow modes from libration to circulation for model systems of varying total mass ratio and initial eccentricities including (Obj)GJ 876(Obj). The attainment of small eccentricities and coupling to the fast mode results in conversion of the slow mode to circulation before the fast mode. Investigations are carried out for stochastic diffusion rates, proportional to the mean square stochastic force amplitude ranging over several orders of magnitude. The life time of the resonant angles is found to be inversely proportional to the diffusion rate except in the case of systems with low total mass in the earth mass range.

In Sect. 5 we exploit the tendency of the slow mode, related to the angle between the apsidal lines, to be driven to circulate while the fast mode still librates in stochastically forced systems, to make a model for the formation of the (Obj)HD 128311(Obj) system which may be in such a state and was not readily understood in terms of convergent migration models for producing the commensurability. We combine the effects of such migration and stochastic forcing, showing that during the migration phase, while the librations tend to be stabilized, the slow mode is readily converted to circulation while the fast mode continues to librate. Good agreement is obtained with the somewhat uncertain observed orbital configuration. Finally in Sect. 6 we summarize and discuss our results.

2 Basic equations

We begin by writing down the equations of motion for a single planet moving in a fixed plane under a general Hamiltonian H in the form (see e.g. Snellgrove et al. 2001; Papaloizou 2003)
 
          $\displaystyle \dot E$ = $\displaystyle -n\frac{\partial H}{\partial \lambda}$ (1)
$\displaystyle \dot G$ = $\displaystyle -\left(\frac{\partial H}{\partial \lambda}+\frac{\partial H}{\partial \varpi}\right)$ (2)
$\displaystyle \dot \lambda$ = $\displaystyle \frac{\partial H}{\partial L} + n \frac{\partial H}{\partial E}$ (3)
$\displaystyle \dot \varpi$ = $\displaystyle \frac{\partial H}{\partial L}\cdot$ (4)

Here the angular momentum of the planet is G and the energy is E. For orbital motion around a central point mass M we have
$\displaystyle G = m\sqrt{ \mathcal{G}Ma(1-e^2)} \hspace{0.5cm} {\rm and}$     (5)
$\displaystyle E = -{{ \mathcal{G}Mm}\over{2a}},$     (6)

where $ \mathcal{G}$ is the gravitational constant, a the semi-major axis and e the eccentricity. The mean longtitude is $\lambda = n (t-t_0) +\varpi ,$ where n the mean motion, with t0 being the time of periastron passage and $\varpi $ being the longitude of periastron.

2.1 Additional forcing of a single planet

In order to study the phenomena such as stochastic forcing we need to consider the effects of an additional external force per unit mass ${\vec F}$ which may not be described using a Hamiltonian formalism. However, as may be seen by considering general coordinate transformations starting from a Cartesian representation, the equations of motion are linear in the components of ${\vec F}$. Because of this we may determine them by considering forces of the form ${\vec F} = (F_x, F_y)$for which the Cartesian components are constant. Having done this we may then suppose that these vary with coordinates and time in an arbitrary manner. Following this procedure we note that when ${\vec F},$ as in the above form is constant, we may derive the equations of motion by replacing the original Hamiltonian with a new Hamiltonian defined through

$\displaystyle H \rightarrow H - m\left( F_x\;x+F_y\;y\right) = H - m\left({\vec r}\cdot{\vec F}\right).$     (7)

The additional terms proportional to the components of ${\vec F}$ correspond to the Gaussian form of the equations of motion (Brouwer & Clemence 1961).

The various derivatives involving ${\vec r}$ can be calculated by elementary means and expressed in terms of $E, G, \lambda$ and $\varpi.$ One thus finds additional contributions to the equations of motion (1)-(4), indicated with a subscript F, in the form

   
                              $\displaystyle \dot G_F$ = $\displaystyle m\left( \frac{\partial }{\partial \lambda } + \frac{\partial }{\p... ...\vec F}\right) = m\;\left( {\vec r}\times {\vec F} \right) \cdot \hat{\vec e}_z$ (8)
$\displaystyle \dot E_F$ = $\displaystyle m n\frac{\partial }{\partial \lambda} \left( {\vec r} \cdot {\vec F} \right) = m\;( {\vec v} \cdot {\vec F})$ (9)
$\displaystyle \dot \varpi_F$ = $\displaystyle - m\frac{\partial }{\partial L} \left( {\vec r} \cdot {\vec F} \right) \ \ \ {\rm or \ \ equivalently}$ (10)
$\displaystyle \dot \varpi_F$ = $\displaystyle \frac{ \sqrt{(1-e^2)} }{n a~e} \left[F_\theta \left( 1 + \frac{1}{{1-e^2}} \frac {r}{a} \right) \sin f - F_r \cos f \right]$ (11)
$\displaystyle \dot \lambda_F$ = $\displaystyle - m\left(\frac{\partial }{\partial L} + n \frac{\partial}{\partial E} \right) \left( {\vec r} \cdot {\vec F} \right)$ (12)
  = $\displaystyle \left( 1- \sqrt{1-e^2}\right) \dot \varpi_F +\frac{2an}{ \mathcal{G}M} \;({\vec r} \cdot {\vec F}) ,$ (13)

where the true anomaly f is defined as the difference between the true longitude and the longitude of periastron, $f = \theta - \varpi.$ Note that from (9) we obtain
$\displaystyle \dot a_F \ \ = \ \ -{2 a \dot n\over 3 n} = \frac {2( {F_r} e \sin f + {F_{\theta}}(1+ e \cos f))} {n\sqrt{1-e^2}},$     (14)

and from (8) together with (9) we obtain
$\displaystyle \dot e_F = \frac{G(2 E \dot G + G \dot E) }{ \mathcal{G}^2m^3M^2 e}\cdot$     (15)

In the limit $e \ll 1$ this becomes (ignoring terms O(e) and smaller)
$\displaystyle \dot e_F = {F_r} \frac 1 {an} \sin f + {F_{\theta}} \frac 1{an} 2 \cos f.$     (16)

Furthermore in this limit we may replace f by $f= \lambda -\varpi = n(t-t_0).$

We remark that the above formalism results in Eq. (11) which gives an expression for $\dot \varpi_F$ that indicates that this quantity diverges for small e as 1/ e. We comment that, as is well known, this aspect results from the choice of coordinates used and is not associated with any actual singularity or instability in the system. This is readily seen if one uses $h=e\sin\varpi,$ and $k=e\cos\varpi$ as dynamical variables rather than e and $\varpi.$ The former set behave like Cartesian coordinates, while the latter set are the corresponding cylindrical polar coordinates. When the former set are used, potentially divergent terms $\propto$1/e do not appear. This can be seen from (11) and (16) which give in the small e limit

 
$\displaystyle \dot h_F = -{F_r} \frac 1 {an} \sin \lambda + {F_{\theta}} \frac 1{an} 2 \cos \lambda$     (17)
$\displaystyle \dot k_F = {F_r} \frac 1 {an} \sin \lambda + {F_{\theta}} \frac 1{an} 2 \cos \lambda .$     (18)

Abrupt changes to $\varpi $ may occur when h and k pass through the origin in the (h,k) plane. But this is clearly just due to a coordinate singularity rather than a problem with the physical system which changes smoothly as this transition occurs. The abrupt changes to the $\varpi $ coordinate occur because very small perturbations to very nearly circular orbits produce large changes to this angle.

2.2 Multiple planets

Up to now we have considered a single planet. However, it is a simple matter to generalize the above formalism so that it applies to a system of two planets. Here we follow closely the discussions in Papaloizou (2003) and Papaloizou & Szuszkiewicz (2005). Excluding stochastic forcing for the time being, we start from the Hamiltonian formalism describing their mutual interactions using Jacobi coordinates (see Sinclair 1975). In this formalism the radius vector ${\vec r}_2,$ of the inner planet of reduced mass m2is measured from Mand that of the outer planet, ${\vec r }_1,$of reduced mass m1 is referred to the centre of mass of M and m2.Thus from now on we consistently adopt a subscripts 1 and 2for coordinates related to the outer and inner planets respectively.

The required Hamiltonian correct to second order in the planetary masses is given by

$\displaystyle H = {1\over 2} ( m_1 \vert \dot {\vec r}_1\vert^2 +m_2\vert \dot ... ...cal{G}m_{1}m_2 {\vec r}_1\cdot {\vec r}_2 \over \vert {\vec r}_{1}\vert^3}\cdot$     (19)

Here M1=M+m1, M2= M + m2 and ${\vec r}_{12}= {\vec r}_{2}- {\vec r}_{1}$. The Hamiltonian can be expressed in terms of $E_i,G_i,\varpi_i, \lambda_i$, i=1,2 and the time t. The energies are given by $E_i = - \mathcal{G}m_iM_{i}/(2a_i),$and the angular momenta $G_i = m_i\sqrt{ \mathcal{G}M_{i}a_i(1-e_i^2)}$ with ai and ei denoting the semi-major axes and eccentricities respectively. The mean motions are $n_i = \sqrt{ \mathcal{G}M_{i}/a_i^3}.$

The Hamiltonian may quite generally be expanded in a Fourier series involving linear combinations of the three angular differences $\lambda_i - \varpi_i, i=1,2$ and $\varpi_1 -\varpi_2$(e.g. Brouwer & Clemence 1961).

Near a first order p+1: p resonance, we expect that both $\phi_1 = (p+1)\lambda_1-p\lambda_2-\varpi_2, $ and $\phi_2 = (p+1)\lambda_1-p\lambda_2-\varpi_1,$will be slowly varying. Following standard practice (see e.g. Papaloizou & Szuszkiewicz 2005; Papaloizou 2003) only terms in the Fourier expansion involving linear combinations of $\phi _1$ and $\phi _2$as argument are retained because only these are expected to lead to large long term perturbations.

The resulting Hamiltonian may then be written in the general form H=E1+E2+ H 12, where

\begin{displaymath}H _{12}= -\frac{ \mathcal{G}m_1m_2}{a_1}\sum C_{k,l}\left( \frac{a_1}{a_2}, e_1,e_2 \right) \cos (l\phi_1 +k\phi_2), \end{displaymath} (20)

where in the above and similar summations below, the sum ranges over all positive and negative integers (k,l) and the dimensionless coefficients Ck,ldepend on e1,e2 and the ratio a1/a2 only. We also replace Mi by M.

2.2.1 Equations of motion

The equations of motion for each planet can now be easily derived that take into account the contributions due to their mutual interactions (see Papaloizou & Szuszkiewicz 2005; Papaloizou 2003) and contributions from (8)-(13). The latter terms arising from external forcing are indicated with a subscript F. We thus obtain to lowest order in the perturbing masses:
     
$\displaystyle {{\rm d}n_1\over {\rm d}t}$ = $\displaystyle {3(p+1)n_1^2 m_2 \over M } \sum C_{k,l}(k+l)\sin (l\phi_1 +k\phi_2)$ $\displaystyle + \left({\rm d}n_1\over {\rm d}t\right)_F$     (21)
$\displaystyle {{\rm d}n_2\over {\rm d}t}$ = $\displaystyle -{3pn_2^2 m_1 a_2\over M a_1 }\sum C_{k,l}(k+l)\sin (l\phi_1 +k\phi_2)$$\displaystyle + \left({\rm d}n_2\over {\rm d}t\right)_F$     (22)
$\displaystyle {{\rm d}e_1\over {\rm d}t}$ = $\displaystyle - {m_2 n_1 \sqrt{1-e_1^2} \over e_1 M } \cdot \sum C_{k,l}\sin (l\phi_1 +k\phi_2)$     (23)
$\displaystyle \cdot \left[k-(p+1)(k+l)\left(1-\sqrt{1-e_1^2}\right)\right] + \left({\rm d}e_1\over {\rm d}t\right)_F$  
$\displaystyle {{\rm d}e_2\over {\rm d}t}$ = $\displaystyle - {m_1 a_2 n_2 \sqrt{1-e_2^2} \over a_1 e_2 M }\cdot \sum C_{k,l}\sin (l\phi_1 +k\phi_2)$  
$\displaystyle \cdot \left[l+p(k+l)\left(1-\sqrt{1-e_2^2}\right)\right] + \left({\rm d}e_2\over {\rm d}t\right)_F$     (24)
$\displaystyle {{\rm d}\phi_2 \over {\rm d}t} = (p+1)n_1- pn_2 - \sum (D_{k,l} +E_{k,l})\cos(l\phi_1 +k\phi_2) + \left({\rm d}\phi_2\over {\rm d}t\right)_F$     (25)
$\displaystyle {{\rm d}\phi_1 \over {\rm d}t} = (p+1)n_1- pn_2 - \sum ( D_{k,l} ... ...k,l} )\cos(l\phi_1 +k\phi_2) + \left({\rm d}\phi_1\over {\rm d}t\right)_F \cdot$     (26)

Here
$\displaystyle D_{k,l} = {2(p+1)n_1a_1^2m_2\over M}{\partial \over \partial a_1}... ...pn_2a_2^2m_1 \over M}{\partial \over \partial a_2}\left( C_{k,l}/a_1 \right)\!,$     (27)
$\displaystyle E_{k,l} = {n_1m_2 \left((p+1)\left(1 -e_1^2\right)-p\sqrt{1-e_1^2... ...t{1-e_2^2}-1+e_2^2\right) \over a_1 e_2 M}{\partial C_{k,l} \over \partial e_2}$     (28)

and
$\displaystyle F_{k,l} = {(p+1)n_1m_2 \left(1 -e_1^2-\sqrt{1-e_1^2}\right) \over... ...-e_2^2\right)\right) \over a_1 e_2 M}{\partial C_{k,l} \over \partial e_2}\cdot$     (29)

Note that $\phi_2 - \phi_1 = \varpi_2 - \varpi_1\equiv \Delta \varpi $is the angle between the two apsidal lines of the two planets. We also comment that, up to now, we have not assumed that the eccentricities are small and that, in additional to stochastic contributions, the external forcing terms may in general contain contributions from very slowly varying disk tides but we shall not consider these further in this article.

2.2.2 Modes of libration

We first consider two planets in resonance with no external forces acting in order to identify the possible libration modes. We then consider the effects of the addition of external stochastic forcing. In the absence of external forces Eqs. (21)-(26) can have a solution for which ai and ei are constants and the angles $\phi _1$ and $\phi _2$ are zero. In general other values for the angles may be possible but such cases do not occur for the numerical examples presented below. When the angles are zero Eqs. (25) and (26) provide a relationship between e1 and e2 (see e.g. Papaloizou 2003).

We go on to consider small amplitude oscillations or librations of the angles about their above equilibrium state. Because two planets are involved there are two modes of oscillation which we find convenient to separate and describe as fast and slow modes. Assuming the planets have comparable masses, the fast mode has libration frequency ${\propto} \!\sqrt{m}$and the slow mode has libration frequency ${\propto} m$. These modes clearly separate as the planet masses are decreased while maintaining fixed eccentricities.

2.2.3 Fast mode

To obtain the fast mode we linearize (21)-(26) and neglect second order terms in the planet masses. This is equivalent to neglecting the variation of Dk,l, Ek,l and Fk,l in Eqs. (25) and (26) which then require that $\phi_1 = \phi_2$ very nearly for this mode. Noting that for linear modes of the type considered here and in the next section, Eqs. (21)-(24) imply that $\dot n_i$ and $\dot e_i$ are proportional to linear combinations of the librating angles, differentiation of either of Eqs. (25) or (26) with respect to time then gives for small amplitude oscillations
$\displaystyle {{\rm d}^2\phi_i \over {\rm d}t^2} +\omega_{lf}^2 \phi_i = 0, \ \ \ ( i= 1,2 ),$     (30)

where

\begin{displaymath}\omega_{lf}^2= - {3p^2n_2^2 m_2 \over M } \left(1+{ a_2 m_1\over a_1m_2}\right)\cdot \sum C_{k,l}(k+l)^2 \end{displaymath} (31)

and we have used the resonance condition that (p+1)n1 = p2 n2which is satisfied to within a correction of order $\sqrt{m_1/M}.$Note that for this mode the fact that $\phi_1 = \phi_2$ very nearly, implies that $\varpi_2 -\varpi_1$ is small. Thus that quantity does not participate in the oscillation.

2.2.4 Slow mode

In this case we look for low frequency librations with frequency  ${\propto} m_1.$ Equations (25) and (26) imply that, for such oscillations, to within a small relative error of order m1/M, (p+1)n1 = pn2 throughout.

Equations (21) and (22) then imply that the two angles are related by $\phi_2=\beta \phi_1,$ where $\beta = - \sum C_{k,l}(k+l)l)/(\sum C_{k,l}(k+l)k).$Subtracting Eq. (26) from Eq. (25), differentiating with respect to time and using this condition results in an equation for $\zeta = \phi_ 2- \phi_1 = \varpi_2 -\varpi_1 = \Delta \varpi$

$\displaystyle {{\rm d}^2\zeta \over {\rm d}t^2} +\omega_{ls}^2 \zeta = 0,$     (32)

where

\begin{displaymath}\omega_{ls}^2= \alpha_1 {n_1m_2\sqrt{1 -e_1^2} \over e_1 M (... ...^2} \over a_1 e_2 M (1-\beta)}{\partial W \over \partial e_2} \end{displaymath} (33)

with
$\displaystyle \alpha_1 = \sum C_{k,l}\left[k-(p+1)(k+l)\left(1-\sqrt{1-e_1^2}\right)\right](k\beta+l),$      
$\displaystyle \alpha_2 = \sum C_{k,l}\left[l+p(k+l)\left(1-\sqrt{1-e_2^2}\right)\right](k\beta+l)$      

and
$\displaystyle % W = \left( {n_1a_2m_2\sqrt{1 -e_2^2} \over a_1 e_2 M}{\partial ... ...1 -e_1^2} \over e_1 M}{\partial \over \partial e_1} \right) \cdot \sum C_{k,l}.$      

Although the expressions for the mode frequencies are complicated, the fast frequency scales as the square root of the planet mass and the slow frequency scales as the planet mass independent of the magnitude of the eccentricities while both scale as the characteristic mean motion of the system. Furthermore although we have considered small amplitude librations and accordingly obtained harmonic oscillator equations, the treatment can be extended to consider finite amplitude oscilations and generalized pendulum equations as long as the two mode frequencies can be separated. However, we shall not consider this aspect further here.

2.2.5 Librations with external forcing

When external forcing is included source terms appear on the right hand sides of Eqs. (30) and (32). We shall assume that the forcing terms are small so that terms involving products of these and both the libration amplitudes and the planet masses may be neglected. Then in the case of the slow mode, repeating the derivation given above including the forcing terms, we find that (32) becomes
$\displaystyle {{\rm d}^2\zeta \over {\rm d}t^2} +\omega_{ls}^2 \zeta = {{\rm d}\over {\rm d}t}\left(\dot \varpi_{2F} - \dot \varpi_{1F} \right)\!.$     (34)

The quantities $\dot\varpi_{iF}$ are readily obtained for each planet from (11). From this we see that for small eccentricities, $\dot\varpi_{iF} \propto 1/e_i,$indicating large effects when ei is small. As already discussed in Sect. 2.1 this feature arises from a coordinate singularity rather than physically significant changes to the system.

A similar description may be found for the fast mode. In this case, neglecting terms of the order of the square of the planet masses or higher, one may use Eqs. (25) and (26) to obtain an equation for $Q \equiv \phi_1 $ in the form

\begin{displaymath}% {{\rm d}^2 Q \over {\rm d}t^2} +\omega_{lf}^2 Q = {{\rm d}... ...eft[ \dot \phi_{1F} \right ] +(p+1)\dot n_{1F}-p\dot n_{2F}. \end{displaymath} (35)

If equations for quantities that are regular functions of the Cartesian like coordinates $(h_i=e_i\sin\varpi_i$ and $k_i=e_i\cos\varpi_i,~ i=1,2)$ are found, the source terms arising from the external forcing are regular as $e_i \rightarrow 0$. To illustrate this we consider the pair
$\displaystyle y\equiv e_1e_2\sin\zeta= h_2k_1-h_1k_2\quad {\rm and}$      
$\displaystyle z\equiv e_1e_2\cos\zeta= k_1k_2-h_1h_2.$      

It is readily seen that their time derivatives are given by
 
$\displaystyle {\dot y}={{\rm d}(e_1e_2)\over {\rm d}t}\sin\zeta + e_1e_2\cos\ze... ...m d}t}\right\vert _U + \left(\dot \varpi_{2F} - \dot \varpi_{1F} \right)\right]$     (36)
$\displaystyle {\dot z}= \cos\zeta\left[\left.{{\rm d}(e_1e_2)\over {\rm d}t}\ri... ...dot e_{1F} + e_1\dot e_{2F}\right]- e_1e_2\sin\zeta{{\rm d}\zeta\over {\rm d}t}$     (37)

here |U denotes evaluation without external forcing. There is now no divergence of the external forcing terms as $e_i \rightarrow 0.$If as above small amplitude slow mode librations of $\zeta $ about zero are considered, given that without external forcing $\dot e_i \propto \zeta,$ the first term on the right hand side of (36) is either quadratic in $\zeta,$ or proportional to the product of the exernal forcing and libration amplitudes and so, as above may be neglected. The same applies to last term on the right hand side of (37). In addition $\cos\zeta$ may be replaced by unity. On taking the time derivatives of (36) and (37), we obtain under the same approximation scheme
 
$\displaystyle {\ddot y} = e_1e_2\left.{{\rm d}^2\zeta\over {\rm d}t^2}\right\ve... ... {\rm d}t}\left[ e_1e_2\left( \dot \varpi_{2F} - \dot \varpi_{1F}\right)\right]$     (38)
$\displaystyle {\ddot z} = e_2\left.{{\rm d}^2 e_1\over {\rm d}t^2}\right\vert _... ... _U +{{\rm d} \over {\rm d}t}\left[ e_1 \dot e_{2F} + e_2\dot e_{1F} \right]\!.$     (39)

Making use of (32) and its counterparts for ${\ddot e_1}$ and ${\ddot e_2}$ for the unforced motion we obtain
 
$\displaystyle {\ddot y} = -\omega_{ls}^2 y +{{\rm d} \over {\rm d}t}\left[ e_1e_2\left( \dot \varpi_{2F} - \dot \varpi_{1F}\right)\right]$     (40)
$\displaystyle {\ddot {\delta{z}}} = -\omega_{ls}^2 \delta z +{{\rm d} \over {\rm d}t}\left[ e_1 \dot e_{2F} + e_2\dot e_{1F}\right]\!,$     (41)

where $\delta z= z-z_0,$ z0 being the value of e1e2 at the centre of the oscillation in z. These may be used as an alternative to (34), and are without potentially singular forcing terms. We note that for the purpose of this paper, both descriptions provide exactly the same physical content.

Equations (35) and (34) (or (40) together with (41)) form a pair of equations for the stochastically forced fast and slow modes respectively. We comment that this mode separation is not precise. However, it can be made so by choosing appropriate linear combinations of the above modes. Numerical results confirm that Q predominantly manifests the fast mode and $\zeta $ or y the slow mode, so we do not expect such a change of basis to significantly affect conclusions.

We further comment that because $\dot\phi_{1F}$ contains $\dot\varpi_{2F}$but not  $\dot\varpi_{1F},$ for small eccentricities there are only potential forcing terms $\propto$1/e2that occur when forcing is applied to the inner planet. As this planet has the larger ecentricity for the situations we consider, small eccentricities are not found to play any significant role in this case.

Each mode responds as a forced oscillator. We suppose the forcing contains a stochastic component which tends to excite the respective oscillator and ultimately convert libration into circulation. But we stress that the above formulation as well as developments below assume small librations, so we may only assess the initial growth of oscillation amplitude. However, inferences based on the structure of the non linear governing equations and an extrapolation of the linear results enable successful comparison with numerical results.

2.3 Stochastic forces

We assume that turbulence causes the external force per unit mass $(F_r, F_{\theta})$ acting on each planet to be stochastic. For simplicity we shall adopt the simplest possible model. Regarding the components of the force, per unit mass expressed in cylindrical coordinates to be a function of time t, we assume any one of these satisfies the relation $F_i(t)F_i(t')= \langle F_i^2 \rangle g(\vert t-t'\vert) $ where the autocorrelation function g(x) is such that $\int^{\infty}_0 g(x) {\rm d}x = \tau_{\rm c},$where $\tau_{\rm c}$ is the correlation time and $\sqrt{\langle F_i^2 \rangle}$ is the root mean square value of the i component. It should be noted that an ensemble average is implied on the left hand side. Again for simplicity we assume that different components acting on the same planet as well as components acting on different planets are uncorrelated. We note that in general the root mean square values as well as $\tau_{\rm c}$ may depend on t, but we shall not take this into account here and simply assume that these quantities are constant and the same for each force component.

We note the stochastic forces make quantities they act on undergo a random walk. Thus if for example $\dot A = F_i$ for some quantity A (note that constants or slowly varying quantities originally multiplying Fi may be absorbed by a redefinition of A and so do not materially affect the discussion given below), the square of the change of A occuring after a time interval t is given by

                         $\displaystyle (\Delta A)^2$ = $\displaystyle \int^t_0 \int^t_0 F_i(t')F_i(t''){\rm d}t'{\rm d}t''$  
  = $\displaystyle \int^t_0 \int^t_0 \langle F_i^2 \rangle g(\vert t'-t''\vert) {\rm d}t'{\rm d}t'' \rightarrow D_i t,$ (42)

here we take the limit where $t/\tau_{\rm c} $ is very large corresponding to an integration time of very many correlation times and $D_i = 2\langle F_i^2 \rangle \tau_{\rm c}$ is the diffusion coefficient.

When the evolution of a stochastically forced planetary orbit is considered, it is more appropriate to consider a model governing equation for A of the generic form

\begin{displaymath}\dot A = F_i \sin(nt), \end{displaymath} (43)

where we recall that $2\pi/n$ is the orbital period (but note that a different value could equally well be considered). We note in passing that, by shifting the origin of time, an arbitrary phase may be added to the argument of the $\sin$ without changing the results given below. One readily finds that Eq. (42) is replaced by
                         $\displaystyle (\Delta A)^2$ = $\displaystyle \int^t_0 \int^t_0 F_i(t')F_i(t'')\sin(n t') \sin(n t'') {\rm d}t'{\rm d}t''$  
  = $\displaystyle \int^t_0 \int^t_0\langle F_i^2 \rangle g(\vert t'-t''\vert) \sin(n t')\sin(n t''){\rm d}t' {\rm d}t''$  
  $\textstyle \rightarrow$ $\displaystyle {\gamma Dt\over 2},$ (44)

where $\gamma(n) = \int^{\infty}_0 g(x)\cos(n x) {\rm d}x .$

Note that when $n\tau_{\rm c}\ll 1,$ corresponding to the correlation time being much less than the orbital period, $\gamma \rightarrow 1.$For larger $\tau_{\rm c},$ $\gamma <1$ gives a reduction factor for the amount of stochastic diffusion. For example if we adopt an exponential form for the autocorrelation function such that

\begin{eqnarray*}g(\vert t'-t''\vert) = \exp \left(- \frac{\vert t'-t''\vert}{\tau_{\rm c}}\right), \end{eqnarray*}


we find
$\displaystyle \gamma = \frac 1{1+n^2 \tau_{\rm c}^2}$     (45)

and for the purposes of comparison with numerical work we shall use this from now on.

2.3.1 Stochastic forcing of an isolated planet

We begin by considering the effect of stochastic forcing on a single isolated planet. In this case we may obtain a statistical estimate for the characteristic growth of the orbital parameters as a function of time by integrating Eqs. (9), (11), (16), (17) and (18) with respect to time directly. We may then apply the formalism leading to the results expressed in generic form through Eqs. (42)-(45) to obtain estimates for the stochastic diffusion of the orbital elements in the limit of small eccentricity in the form
    
$\displaystyle (\Delta a)^2 = 4 \frac { D t }{ n^2 }$     (46)
$\displaystyle (\Delta e )^{2} = 2.5 \frac {\gamma D t }{ n^2 a^2 }$     (47)
$\displaystyle (\Delta \varpi)^2 = \frac{2.5}{e^2} \frac{\gamma Dt}{n^2 a^2}$     (48)
$\displaystyle (\Delta h )^{2} = 2.5 \frac {\gamma D t }{ n^2 a^2 }$     (49)
$\displaystyle (\Delta k )^{2} = 2.5 \frac {\gamma D t }{ n^2 a^2 }\cdot$     (50)

We note again that the 1/e2 dependence of $(\Delta \varpi)^2$which arises from the coordinate singularity discussed in Sect. 2.1 does not appear for $(\Delta h )^{2}$ and $ (\Delta k )^{2}.$Note that the definition of (h,k) imply consistently with the above that
 
$\displaystyle (\Delta h )^{2} = (\Delta e )^{2}\sin^2\varpi+ e^{2} (\Delta \varpi)^2\cos^2\varpi$     (51)
$\displaystyle (\Delta k )^{2} = (\Delta e )^{2}\cos^2\varpi+ e^{2} (\Delta \varpi)^2\sin^2\varpi.$     (52)

  2.3.2 Stochastic variation of the resonant angles in the two planet case

We now consider the effects of stochastic forcing on the resonant angles. The expressions are more complicated than in the previous section because there are more variables involved. However, we basically follow the formalism outlined in Sect. 2.3.

While the libration amplitude is small enough for linearization to be reasonable, the evolution is described by Eqs. (35) and (34) (or (40) together with (41)). These may be solved by the method of variation of parameters. Assuming the amplitude is zero at t=0, the solution of Eq. (40) is given by

$\displaystyle y = \sin(\omega_{ls}t)\int^t_0{S\!_y\cos(\omega_{ls}t)\over \omeg... ...- \cos(\omega_{ls}t)\int^t_0{S\!_y\sin(\omega_{ls}t)\over \omega_{ls}}{\rm d}t,$     (53)

where $S\!_y= {\rm d}(e_1e_2(\dot\varpi_{2F}- \dot \varpi_{1F}))/{\rm d}t.$ There are corresponding expressions that can be obtained from Eqs. (41) and (35) for $\delta z$ and Q respectively.

Equation (53) may be regarded as describing a harmonic oscillator whose amplitude varies in time such that the square of the amplitude after a time interval t is given by

\begin{displaymath}(\Delta y )^2=\left(\int^t_0{S\!_y\sin(\omega_{ls}t)\over \om... ..._0{S\!_y\cos(\omega_{ls}t)\over \omega_{ls}}{\rm d}t\right)^2. \end{displaymath} (54)

The corresponding expression from Eq. (41) is

\begin{displaymath}(\Delta \delta z )^2=\left(\int^t_0{S\!_z\sin(\omega_{ls}t)\o... ...0{S\!_z\cos(\omega_{ls}t)\over \omega_{ls}}{\rm d}t\right)^2, \end{displaymath} (55)

where $S\!_z = {\rm d}(e_1\dot e_{2F}+e_1\dot e_{2F})/{\rm d}t$. And finally the equation obtained from Eq. (35) is

\begin{displaymath}(\Delta Q)^2 = \left(\int^t_0{S\!_Q\sin(\omega_{lf}t)\over \o... ..._0{S\!_Q\cos(\omega_{lf}t)\over \omega_{lf}}{\rm d}t\right)^2, \end{displaymath} (56)

where $S\!_Q= {\rm d}(\dot \phi_{1F})/{\rm d}t +(p+1)\dot n_{1F} - p\dot n_{2F}.$

We now evaluate the expectation values of these using the formalism of Sect. 2.3. For simplicity and as considered numerically later we shall specialize to the case when stochastic forces act only on the outer planet (but see Sect. 6 below). Taking Eq. (54), we perform an integration by parts neglecting the end point contribution as these are associated with subdominant contributions increasing less rapidly than t for large t, to obtain

                   $\displaystyle (\Delta y )^2$ = $\displaystyle (\Delta (e_1e_2)\sin\zeta )^2$  
  = $\displaystyle \left\vert\int^t_0{e_1e_2\dot \varpi_{1F}\sin(\omega_{ls}t)}{\rm d}t\right\vert^2$  
    $\displaystyle + \left\vert\int^t_0{e_1e_2\dot\varpi_{1F}\cos(\omega_{ls}t)}{\rm d}t\right\vert^2.$ (57)

We then find, working in the limit of a small eccentricity e1 (but not necessesarily e2), from the corresponding equation for $\delta z$ that $(\Delta \delta z )^2 = (\Delta y )^2 .$

In dealing with Eq. (56) we neglect $\dot\phi_{1F}$ in S because after integration by parts this leads to a contribution on the order $\omega_{lf}/n$ smaller than that derived from $\dot n_{1F}.$ Thus we simply obtain

$\displaystyle {(\Delta Q)^2\over (p+1)^2} = \left(\int^t_0{\dot n_{1F}\sin(\ome... ...left(\int^t_0{\dot n_{1F}\cos(\omega_{lf}t)\over \omega_{lf}}{\rm d}t\right)^2.$     (58)

We now follow the procedures outlined in Sect. 2.3 obtaining
 
       $\displaystyle (\Delta e_1e_2\sin\zeta )^2$ = $\displaystyle 2.5 {De_2^2 \gamma_s t\over 2a_1^2 n_1^2 }\quad \quad {\rm and}$ (59)
$\displaystyle {(\Delta Q)^2\over (p+1)^2}$ = $\displaystyle {9 D \gamma_ft\over a_1^2\omega_{lf}^2}$ (60)

where
                          $\displaystyle \gamma_f$ = $\displaystyle \frac1{1+\omega_{lf}^2\tau_{\rm c}^2} \quad \quad {\rm and}$ (61)
$\displaystyle \gamma_s$ = $\displaystyle \frac1{1+(n_1+\omega_{ls})^2\tau_{\rm c}^2} +\frac1{1+(n_1-\omega_{ls})^2\tau_{\rm c}^2 }\cdot$ (62)

2.3.3 Growth of libration amplitudes

Equations (59) and (60) express the expected growth of the resonant angle libration amplitudes as a function of time. We remark that these expressions can be simply related to those obtained for a single planet. Thus Eqs. (46) and (60) applied to the outer planet imply that

$\displaystyle (\Delta Q)^2 /(\Delta a_1)^2= 9(p+1)^2 n_1^2\gamma_f /(4 a_1^2\omega_{lf}^2).$      

As we are interested in the case p=1, the width of the libration zone is ${\sim} a_1\omega_{lf}/n_1,$ we see that the time for $(\Delta Q)^2$ to reach unity is comparable for the semi-major axis to diffuse through the libration zone.

Similarly, for small amplitude librations about fixed e2 Eq. (59) gives almost the same result for $(\Delta e_1\sin\zeta )^2$ or $e_2^{-2}(\Delta \delta z)^2$ as that obtained for $(\Delta h )^{2}$ or $(\Delta k)^2$ from Eq. (49) applied to the outer planet. This indicates that $\zeta $, being the angle between the apsidal lines of the two planets, diffuses in the same way as for an isolated outer planet subject to stochastic forces. Thus, in the small amplitude regime, the way this diffusion occurs would appear to be essentially independent of whether the planets are in resonance (but see below).

An important consequence of Eq. (59) is the behaviour of $\zeta $ for small e1. The latter quantity was assumed constant in the analysis. As implied by the discussion of Sect. 2.1 abrupt changes to $\zeta $ are expected when (y,z) passes through the origin. Then even an initially small amplitude libration is converted to circulation. Thus if e1 is small then Eq. (59) indicates that a time $ t \sim 4a_1^2 n_1^2 e_1^2 / (5 D \gamma_s \zeta^2 )$, is required to convert libration to circulation. This can be small if e1 is small. Even if e1 is not small initially, it is important to note that it also undergoes stochastic diffusion (see Eq. (47)) as well as oscillations through its participation in libration. Should e12attain very small values through this process, then from (59) we expect the onset of a rapid evolution of $\zeta.$ Accordingly the attainment of circulation for this angle, should be related to the diffusion of e12 allowing very small values of that quantity to be attained, rather than the direct excitation of libration amplitude. This is particularly the case when e12 starts from relatively small values.

In fact, application of (59) and (60) to the numerical examples discussed below adopting the initial orbital elements indicates that the diffusion of $\zeta $ is significantly smaller than Q unless e1 starts out with a very small value. This would suggest that Q reaches circulation before $\zeta.$ However, this neglects the coupling between the angles that occurs once the libration amplitudes become significant. It is readily seen that it is not expected that Q could circulate while $\zeta $ remains librating as it was initially. One expects to recover standard secular dynamics for $\zeta $ from the governing Eqs. (21)-(25), when these are averaged over an assumed $Q=\phi_1$ circulating with constant $\dot Q.$ As a libration of the initial form would not occur under those conditions, we expect, and find, large excursions or increases in the libration amplitude of $\zeta $ to be correlated with increases in the libration amplitude of Q. This in turn increases the oscillation amplitude of the eccentricity e1, allowing it to approach zero. The consequent rapid evolution of $\zeta $then enables it to pass to circulation. Thus the breaking of resonance is ultimately found to be controlled by the excitation of large amplitude librations for $Q =\phi_1,$which induce $\zeta $ to pass to circulation somewhat before Q itself.

3 Numerical simulations

We have performed numerical simulations of one and two planet systems that allow for the incorporation of additional stochastic forces with the properties described above. These in turn provide a simple prescription for estimating the effects of stochastic gravitational forces produced by density fluctuations associated with disk turbulence. Results have been obtained using both a fifth order Runge Kutta method and the Bulirsch Stoer method (Press et al. 1992; Stoer & Bulirsch 2002) with fixed as well as adaptive timesteps. We have checked that results are converged and thus do not actually depend on the integrator used.

First we discuss the expected scaling of the stochastic forces with the physical parameters of the disk and their implementation in the n-body integrations. In order to clarify the physical mechanisms involved and to check the analytic predictions for stochastic diffusion given by Eqs. (46)-(49) we consider simulations of a single planet undergoing stochastic forcing first. We then move on to consider two planet commensurable systems with and without stochastic forcing. We focus on the way a 2:1 commensurability, corresponding to p=1 is disrupted, highlighting the various evolutionary stages a system goes through as it evolves from a state with a strong commensurability affecting the interaction dynamics, to one where the commensurability is completely disrupted and in some cases a strong scattering occurs. We consider a range of different planet masses and eccentricities (see Table 1).

3.1 Stochastic forces

In order to mimic the effects of turbulence, for example produced by the MRI, it is necessary to calibrate these forces with reference to MHD simulations. As described above, the basic parameters characterising the prescription for stochastic forcing that we have implemented are the root mean square value of the force components per unit mass (in cylindrical coordinates) $\sqrt{\langle F_i^2 \rangle}$ and the auto correlation time $\tau_{\rm c}$.

Table 1:   Parameters of the model systems considered.

From our analytic considerations, we concluded the stochastic forces make the orbital parameters undergo a random walk that is dependent on the force model primarily through the diffusion coefficient $D = 2\langle F_i^2 \rangle \tau_{\rm c} .$

For planets under the gravitational influence of a protoplanetary disk, the natural scale for the force per unit mass components, Fi, is $F_0(r) = \pi \mathcal{G}\Sigma(r)/2$, where $\Sigma$ is the characteristic disk surface density (see e.g. Papaloizou & Terquem 2006). We comment that F0(r) is the gravitational force per unit mass due to a small circular disk patch of radius $r_{\Sigma}$ at a distance $\sqrt{2} \; r_{\Sigma}$ from its centre assuming that all its mass is concentrated there. The result is independent of $r_{\Sigma}.$ The natural correlation time $\tau_{\rm c}$ is the inverse of the orbital angular frequency $\tau_{\rm c,0} = \Omega^{-1}$. To set the natural scale for D, we adopt a minimum mass solar nebula model (MMSN, see Weidenschilling 1977) with

$\displaystyle \Sigma(r) = 4200 \frac{{\rm g}}{{\rm cm}^2} \left(\frac{r}{{\rm 1~ AU}}\right)^{-3/2}\cdot$     (63)

This provides a natural scale for D as a function of the local disk radius and the central stellar mass through
$\displaystyle D_0 = 2 C F_0^2 \tau_{\rm c,0} \ \ = 25 \; C \; \frac{ {\rm cm}^2... ...{1~{\rm AU}}\right)^{-\frac32} \left(\frac{M_*}{1 ~M_\odot}\right)^{-\frac 12},$     (64)

where C is a dimensionless constant.

There are several very uncertain factors which contribute to determining an appropriate value of the dimensionless constant C: the density fluctuations found in MRI simulations are typically ${\delta \rho}/{\rho} = {\delta \Sigma}/{\Sigma} \approx 0.1$ (e.g. Nelson 2005). The presence of a dead zone in the mid plane regions of the disk, where the MRI is not active, has been found to cause reductions in the magnitude of F0 by one order of magnitude or more, as compared to active cases (Oishi et al. 2007). Massive planets open a gap in the disk. Oishi et al. (2007) found that most of the contribution to the stochastic force comes from density fluctuations within a distance of one scaleheight from the planet. When a gap forms, this region is cleared of material leading to the expectation of a substantial decrease in the the magnitude of turbulent density fluctuations. Consequently F0 should be reduced on account of a lower ambient surface density. A factor of $\frac 1 {10}$ seems reasonable although it might be even smaller (de Val-Borro et al. 2006). The correlation time $\tau_{\rm c}$ is actually found to be approximately  $0.5\Omega^{-1}$ (Oishi et al. 2007; Nelson & Papaloizou 2004).

If it is appropriate to include reduction factors to account for all of the above effects, one finds $C= 5\times 10^{-7}$ and we expect a natural scale for the diffusion coefficient to be specified through

\begin{displaymath}D_0 \rightarrow 10^{-5}\left( \frac{r}{1~ \rm AU}\right)^{-3/... ...\odot}\right)^{-1/2}~ \frac{ {\rm cm}^2}{{\rm s}^{3}}\! \cdot \end{displaymath} (65)

The same value of D may be equivalently scaled to the orbital parameters of the planets without reference to the disk by writing
$\displaystyle D_0 = 3.5 \times 10^{6}\left( \frac{r}{1~{\rm AU}}\right)^{-5/2}\... ...ngle \Omega \tau_{\rm c}}{(GM_*)^2}\right) \frac{ {\rm cm}^2}{{\rm s}^{3}}\cdot$     (66)

Thus a value D0=10-5 in cgs units corresponds to a ratio of the root mean square stochastic force component to that due to the central star of about ${\sim} 10^{-6}$ for a central solar mass at 1 AU. It is a simple matter to scale to other locations.

Of course we emphasize that the value of this quantity is very uncertain, a situation that is exacerbated by its proportionality to the square of the magnitude of the stochastic force per unit mass. For this reason we perform simulations for a range of Dcovering many orders of magnitude.

3.2 Numerical implementation

The procedure we implemented, uses a discrete first order Markov process to generate a correlated noise that is continuous and added as an additional force. The Markov process is a statistical process which is defined by two parameters, the root mean square of the amplitude and the correlation time $\tau_{\rm c}$ (Kasdin 1995). It has a zero mean value and has no memory. This has the advantage that previous values do not need to be stored. The autocorrelation function decays exponentially and thus mimics the autocorrelation function measured in MHD simulations by Oishi et al. (2007).

3.3 Stochastic forces acting on a single planet

We first investigate the long term effect of stochastic forces on a single isolated planet. The initial orbital parameters were taken to be the the observed parameters of (Obj)GJ 876 b(Obj) (see Table 1) and the central star had a mass of $0.38~{M}_{\odot}$. In this simulation, we use the reduced diffusion coefficients
                          D = $\displaystyle 8.2 \times 10^{-5} ~\frac{ {\rm cm}^2}{{\rm s}^{3}} \quad {\rm and}$ (67)
D = $\displaystyle 8.2 \times 10^{-3} ~\frac{ {\rm cm}^2}{{\rm s}^{3}}$ (68)

which can be represented by a correlation time of half the orbital period and a specific force with root mean square values $\sqrt{\langle F^2\rangle }=4.05\times10^{-6}~{\rm cm}/{\rm s}^{2}$ and $\sqrt{\langle F^2\rangle }=4.05\times10^{-5}~{\rm cm}/{\rm s}^{2}$, respectively (see also Eq. (66) above).

The resulting random walks undergone by e, a and $\varpi $ are plotted in Fig. 1 for six different realisations for each of the two diffusion parameters. The spreading rates can be estimated from Eqs. (46)-(49). These analytic predictions are plotted as solid lines. Clearly the numerical model is in broad agreement with the random walk description with a spreading that scales with D as expected. We have performed simulations for a variety of D and get results that are fully consistent with this in all cases.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330f1a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f1b.eps}\end{figure} Figure 1:

Time evolution of the longitude of periastron $\varpi $, the eccentricity e and the semi major axis a for a single planet. The initial orbital parameters of the planet were taken to be those of (Obj)GJ 876 b(Obj) (see Table 1). Six different realisations starting from the same initial conditions are shown in each panel. The diffusion coefficients are $D = 8.2 \times 10^{-5} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the upper three panels and $D = 8.2 \times 10^{-3} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the lower three respectively. The solid lines correspond to the analytic predictions for the amount of spreading (see text).

Open with DEXTER

3.4 Illustration of the modes of libration in a two planet system

We now go on to consider two planet systems. As an illustrative example we consider the (Obj)GJ 876(Obj) system (see Table 1). Note that we can easily scale all physical quantities and extend the discussion to other systems (see Sect. 6 for a discussion). We begin by considering the evolution of the system without stochastic forcing in order to characterize the modes of libration of the resonant angles and other orbital parameters. In particular we identify the fast and slow modes discussed in Sect. 2.2.2.

The time evolution of the resonant angles and the eccentricities is plotted in Fig. 2. Clearly visible are the slow and fast oscillation modes. The fast mode, which is seen to have a period of about 1.4 years, dominates the librations of e2, and $\phi _1$ while also being present in those of $\phi_2.$ On the other hand the slow mode, which is seen to have a period of about 10 years, dominates the librations of $\zeta $ while also being present in those of e1 and $\phi_2.$

We emphasize the fact that the eccentricities of the two planets participate in the librations and so are not constant. In particular, the eccentricity of (Obj)GJ 876 b(Obj) oscillates around a mean value of 0.03 with an amplitude $\Delta e \approx 0.01 {-} 0.02$. This behaviour involving the attainment of smaller values of the eccentricity has important consequences for stochastic evolution as discussed above (see Sect. 2.3.3) and see also below.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig2.eps} \par\end{figure} Figure 2:

Time evolution of the resonant angles and the eccentricities in the system (Obj)GJ 876(Obj) without turbulent forcing. The dominance of the fast mode with period ${\sim } 1.4$ years in the oscillations of $\phi _1$, and the dominance of the slow mode with period ${\sim } 10$ years in $\zeta $ can be clearly seen.

Open with DEXTER

We remark that similar behaviour occurs for all the systems we have studied, these have a wide range of eccentricities and planet masses. Indeed we note that the period scalings of the fast and slow mode periods with the planet masses are provided by Eqs. (31) and (33) respectively. These indicate, as confirmed by our simulations, that if both masses are reduced by a factor $\Lambda$then the period of the fast mode scales as $\sqrt{\Lambda}$and the period of the slow mode scales as $\Lambda.$

4 Two planets with stochastic forcing

We now consider systems of two planets with stochastic forcing. For simplicity we begin by applying forcing only to the outer planet. We have found that adding the same form of forcing to the inner planet tends to speed up the evolution by approximately a factor of two without changing qualitative details. For illustrative purposes we again start with the (Obj)GJ 876(Obj) system and consider two diffusion coefficients $D = 8.2 \times 10^{-5} ~{{\rm cm}^2}/{{\rm s}^{3}}$ and $D = 8.2\times10^{-4} ~ {\rm cm}^2/{\rm s}^{3}$. We remark that all of our simulations were run with constant values of $D=2\langle F_i\rangle ^2\tau_{\rm c}.$This was done maintaining the parameters $\langle F_i\rangle $ and $\tau_{\rm c}$ to be constant, with $\tau_{\rm c}$ being determined for the initial location of the outer planet. For the cases considered here, there is little orbital migration so this is not a significant feature. Note also that different realizations of the system are likely to become unstable and scatter for higher values of D (see below).

The time evolution of the eccentricities is plotted in Fig. 3. We see fast oscillations superimposed on a random walk. The amplitude of the oscilations, as well as the mean value, changes with time.

Our simple analytic model, assumes slowly changing background eccentricities and semi-major axes and accordingly does not incorporate the oscillations of the eccentricity due to the resonant interaction of the planets. In order to make a comparison, we perform a time average over many periods to get smoothed quantities whose behaviour we can compare with that expected from Eqs. (46)-(49). When this procedure is followed, the evolution is in reasonable accord with that expected from the analytic model provided that allowance is made for the importance of small values of e1 in determining the growth of the libration amplitude of the angle between the apsidal lines of the two planets (see Eq. (49)). The presence of this feature results in the behaviour of the libration amplitude being more complex than that implied by a process governed by a simple diffusive random walk.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330f3a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f3b.eps}\end{figure} Figure 3:

Time evolution of the eccentricities in the (Obj)GJ 876(Obj) system with turbulent forcing included. The diffusion coefficient is $D = 8.2 \times 10^{-5} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the two upper plots and $D = 8.2\times10^{-4} ~ {\rm cm}^2/{\rm s}^{3}$ for the two lower plots. The first (uppermost) and the third plots each show a single run. The second and fourth (lowermost) plots show the time averaged eccentricity for four runs. The averaging interval is 1000 years.

Open with DEXTER

4.1 Disruption of a resonance in stages

Systems with mean motion commensurabilities can be in many different configurations. Here we describe the important evolutionary stages as they appear in a stochastically forced system that starts from a situation in which all the resonant angles show small amplitude libration. For example observations of (Obj)GJ 876(Obj) suggest that the system is currently in such a state with the ratio of the orbital periods P1/P2 oscillating about a mean value of 2.

4.1.1 Attainment of circulation of the angle between the apsidal lines

When the initial eccentricity of the outer planet, e1,is small the excitation of the libration amplitudes of the resonant angles readily brings about a situation where e1 attains very small values. This causes the periastron difference $\zeta $ to undergo large oscillations and eventually circulation (see Eq. (49)). Should stochastic forces cause e1 to reach zero, on account of the coordinate singularity $\varpi_1$ becomes undefined. Subsequently very small perturbations are able to produce a small eccentricity with $\zeta $ undergoing large amplitude librations or circulation (see below). But note that as also mentioned in Sect. 2.1, this occurs without a large physical perturbation to the system. Thus the occurence of this phenomenon does not imply the system ceases to be in a commensurability.

In this context we note that the eccentricity of (Obj)GJ 876 b(Obj) initially is such that $e_1 \sim 0.03$ with values $e_1\sim 0.01$ often being attained during libration cycles. Thus only a small change may cause the above situation to occur. In all cases we have considered, we find that $\zeta $ enters circulation prior to the fast angle $\phi_1,$ which may remain librating until that too is driven into circulation.

4.1.2 Attainment of circulation of the fast angle

Both before and after $\zeta $ enters circulation, stochastic forcing acts to increase the libration amplitude of the fast mode. This mode dominates both the librations of the resonant angle $\phi _1$ and the semi major axes. Eventually $\phi _1$ starts circulation. Shortly afterwards commensurability is lost and P1/P2 starts to undergo a random walk with a centre that drifts away from 2. Note that it is possible for some realizations to re-enter commensurability. For systems with the masses of the observed (Obj)GJ 876(Obj) system, the most likely outcome is a scattering event that causes complete disruption of the system.

4.1.3 A numerical illustration

In order to illustrate the evolutionary sequence described above we plot results for two realisations of the evolution of the (Obj)GJ 876(Obj) system in Fig. 4. For these runs we adopted the diffusion coefficient $D = 0.42~ {\rm cm}^2/{\rm s}^{3}.$ In this context we note that reducing D increases the evolutionary time which has been found, both analyticaly and numerically to be ${\propto}1/D$ (see below).

The times at which the transition from libration to circulation occurs for both the slow and fast angles are indicated by the vertical lines in Fig. 4.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330f4a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f4b.eps}\end{figure} Figure 4:

Time evolution of the resonant angles, the period ratio P2/P1 and the eccentricity, e1, in the (Obj)GJ 876(Obj) system with stochastic forcing corresponding to $D = 0.42~ {\rm cm}^2/{\rm s}^{3}$. The vertical lines indicate when the angles enter circulation for a prolonged period. The realisation illustrated in the lower panel scatters shortly after $\phi _1$ goes into circularization.

Open with DEXTER
Several of the features discussed above and in Sect. 2.3.3 can be seen in Fig. 4. In particular the tendency for the occurence of very small values of e1 to be associated with transitions to circulation of $\zeta $ can be seen for the realisation plotted in the lower panels at around t=80 years and t= 500 years. Such episodes always seem to occur when the libration amplitude of the fast angle $\phi _1$ is relatively boosted, indicating that this plays a role in boosting the slow angle. If the period of time for which e1 attains small values is small and $\phi _1$ recovers small amplitude librations, the slow angle returns to circulation. Thus the attainment of long period circulation for the slow angle is related to the diffusion time for $\phi_1.$We also see from Fig. 4 that the angle $\phi_2,$ which has a large contribution from the slow mode, behaves in the same way as $\zeta $as far as libration/circularization is concerned. We have verified by considering the results from the simulations of GJ 876 LM HE, which started with a larger value of e1, that, as expected, the attainment of circulation of the slow angle takes relatively longer in this case, the time approaching more closely the time when $\phi _1$attains circulation. Also as expected, the time when $\phi _1$ attains circulation is not affected by the change in e1.

4.2 Dependence on the diffusion coefficient

We now consider the stability of the systems listed in Table 1 as a function of D.These systems have a variety of masses and orbital eccentricities. In particular, in view of the complex interaction between the resonant angles discussed above we wish to investigate whether the mean amplitude growth at a given time is indeed proportional to D. As also mentioned above, the value of the diffusion parameter D,that should be adopted, is very uncertain. We have therefore considered values of D ranging over five orders of magnitude. However, the correlation time $\tau_{\rm c}$ is always taken to be given by $\tau_{\rm c} = 0.5\Omega^{-1}$ while the rms value of the force is changed. In order to determine the ``lifetime'' of a resonant angle, we monitor whether it is librating or circulating. Numerically, libration is defined to cease at the first time the angle is seen to reach absolute values largen than 2. We note that the angle can in general regain small values afterwards. However, this is a transient effect and changes the lifetime by no more than a factor of $\sim $2 in all our simulations.

In this context we consider the fast angle $\phi _1$ and the slow angle $\zeta,$though, as we saw above, the latter can be replaced by $\phi _2$which exhibits the same behaviour. As $\phi _1$ is the last to start circulating the resonance is defined to be broken at that point.

Equations (60) and (59) estimate the spreading of the resonant angles as a function of time. We determine the times to attain circulation as the times to attain $(\Delta \phi_i)^2 = 4, i = 1,2 $assuming the initial values for the orbital elements. We plot both the numerical and analytical results in Fig. 5. To remove statistical fluctuations and obtain a mean spreading time, the numerical values for a particular value of D were obtained by averaging over 60 realizations.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330f5a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f5b.eps} \end{figure} Figure 5:

Average time until circulation of the resonant angles $\phi _1$ ( top) and $\phi _2$ ( bottom) in the GJ 876 (indicated with +) , GJ 876LM (indicated with $\times $), GJ 876SE (indicated with $\star $) and GJ 876E (indicated with $\Box $) systems as a function of the stochastic forcing diffusion coefficient D.In the case of the upper panel, the analytic curves explained in the text, are from top to bottom for the GJ 876, GJ 876LM, GJ 876SE and GJ 876E systems respectively. In the lower panel the analytic curves, as explained in the text, are from top to bottom for the GJ 876 (solid curve and top dashed curve), GJ 876LM, GJ 876SE and GJ 876E systems respectively.

Open with DEXTER

From Fig. 5 it is apparent that the evolutionary times scale is ${\propto}1/D$for D varying by many orders of magnitude. The analytic estimates for the libration survival times of $\phi_1,$dominated by the fast mode, obtained using Eq. (60) adopting the initial orbital elements and with the fast libration frequency determined from the simulations, are plotted in the upper panel of Fig. 5. These are in good agreement with the numerical results. However, the estimate for $\phi _2$ based on Eq. (59) using the initial value of e1 overestimates the lifetime by a factor of at least $\sim $40 (see solid line in the lower plot of Fig. 5). Futhermore (59) has no dependence on planet mass which is in clear conflict with the numerical results. This, as discussed above is due to the temporary attainment of small eccentricities. This causes the disruption of the libration of $\zeta $ earlier than would be predicted assuming e1 is constant. In fact this disruption occurs at times that can be up to a factor of 10 times shorter than those required to disrupt $\phi _1$. We estimate the lifetime of librations of $\phi _2$ by calculating the time $\phi _1$ needs to reach values close to 1 (see dashed lines in the lower plot of Fig. 5). As explained above in Sect. 2.3.3, large amplitude variations of $\phi _1$are expected to couple to and excite the slow mode. Variations induced in the eccentricity e1 allow this to reach zero and we lose libration of $\phi_2.$ These simple estimates are in good agreement with the numerical results and accordingly support the idea that $\phi _1$ and $\phi _2$ are indeed coupled in the non linear regime. For low mass planets this simple analytic prediction underestimates the lifetime of $\phi _2$ by a factor of $\sim $2, suggesting that the coupling between the two modes depends slightly on the planet masses.

To confirm our understanding of these processes, we have repeated the calculations with systems that are in the same state as the system discussed above but have higher eccentricities (see GJ 876 LM HE and GJ 876 SE HE in Table 1). The lifetime of the librating state of $\phi _1$ is unchanged, as this does not depend significantly on the eccentricity. However, the lifetime of the librating state of $\phi _2$ is a factor of 2-3 longer. This is due to the fact that a larger excitation of $\phi _2$ is needed to make e1 reach small values in these simulations.

The random walk description breaks down completely when the anticipated disruption time becomes shorter than the libration period. This is expected to happen for small enough masses as can be verified from Fig. 5. It is because the disruption time decreases linearly with the planet masses while the libration period increases as the square root of the planet masses. Then we cannot average over many libration periods. This situation is apparent in Fig. 5 for very short disruption times of the order 100-1000 years where survival times cease to vary linearly with 1/D.

5 Formation of HD 128311

We here discuss the application of the ideas discussed above to understanding the orbital configuration of the planetary system (Obj)HD 128311(Obj). This is (with 99% confidence) in a 2:1 mean motion resonance with the angle $\phi _1$ librating and the angle $\zeta $ circulating so that there is no apsidal corotation (Vogt et al. 2005). However, the orbital parameters are not well constrained. The orgininal Keplerian fit by Vogt et al. (2005) is such that the planets undergo a close encounter after only 2000 years. The values in Table 1 have been obtained from a fit to the observational data that includes interactions between the planets. The values quoted have large error bars. For example the eccentricities e1 and e2 have an uncertainty of 33% and 21%, respectively. Although the best fit doesn't manifest apsidal corotation, the system could undergo large amplitude librations and still be stable.

According to Lee & Peale (2002) the planets should have apsidal corotation if the commensurability was formed by the two planets undergoing sufficently slow convergent inward migration, and if they then both migrated inwards significantly while maintaining the commensurability (see also Snellgrove et al. 2001). Accordingly Sándor & Kley (2006) suggested a possible formation scenario with inward migration as described above, but with an additional perturbing event, such as a close encounter with an additional planet occuring after the halting of the inward migration. This perturbation is needed to alter the behaviour of $\zeta,$ so that it undergoes circulation rather than libration and thus producing orbital parameters similar to the observed ones.

We showed above that stochastic forcing possibly resulting from turbulence driven by the MRI readily produces systems with commensurabilities without apsidal corotation if the eccentricities are not too large. This suggests that a scenario that froms the commensurability through disk induced inward convergent migration might readily produce commensurable systems without apsidal corotation if stochastic forcing is included. Such scenarios are investigated in this section.

We present simulations of the formation of (Obj)HD 128311(Obj) that include migration and stochastic forcing due to turbulence but do not invoke special perturbation events involving additional planets. We find that model systems with orbital parameters resembling the observed ones are readily produced that give better matches than so far provided by Sándor & Kley (2006). The planets in this system are of the order of several Jupiter masses and the eccentricity of one planet can get very small during one libration period. The observed orbital parameters are given in Table 1 and plotted on the right hand side of Figs. 6-8. The mass of the star is $0.8~M_{\odot}$.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig6.eps} \end{figure} Figure 6:

The plots on the left hand side show a possible formation scenario of (Obj)HD 128311(Obj). We plot the observed system on the right hand side as a comparison (see Table 1 and text). The plots show the resonant angles $\phi _1, \phi _2, \zeta $, the eccentricities e1, e2 and period ratio P1/P2 for formation scenarios including turbulence and migration. Resonance capturing occurs after 2000 years.

Open with DEXTER

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig7.eps} \end{figure} Figure 7:

The plots on the left hand side show another possible formation scenario of (Obj)HD 128311(Obj). Again, we plot the observed system on the right hand side as a comparison (see Table 1 and text). The migration timescales in this run are $\tau_{a,1}=16~000$, $\tau_{a,2}=-40~000$ and K=5.5. Note that the eccentricities are larger compared to Fig. 6 because K is smaller.

Open with DEXTER

5.1 Migration forces

We incorporate the non-conservative forces excerted by a protostellar disk that lead to inward migration by following the approach of Lee & Peale (2001) and Snellgrove et al. (2001). In this the migration process is characterized by migration and eccentricity damping timescales, $\tau_a = {a}/{\dot a}$ and $\tau_e = {e}/{\dot e}$. We also define the ratio of the above timescales to be $K=\left\vert\tau_a\right\vert/\tau_e$. This ratio determines the eccentricities in the state of self similar inward migration of the commensurable system that is attained after large times. We verify this result for two different values of K (see below). The procedure is now widely used (e.g. Sándor & Kley 2006; Crida et al. 2008) and we have checked our code against the results of their work.

In this work we allow the planets to form a commensurability through convergent inward migration but stochastic forcing is included, the disk is then removed through having its surface density reduced to zero on a 2000 year timescale. This simultaneously reduces both the migration forces and the stochastic forcing to zero. This procedure is very similar to that adopted by the above authors but we have included stochastic forcing and removed the disk on longer time scales.

The stochastic forces have to have the right balance with respect to the migration rate. We have found that inward migration imparts stability to the resonant system. If the migration rate is too fast relative to the stochastic forcing the migration keeps down the libration amplitudes and we do not get circulation. On the other hand large eccentricity damping favours broken apsidal corotation. This might look counterintuitive at first sight, but remember that the diffusion of $\zeta $ depends on 1/e21,2. Due to the stochastic nature of the problem, it is hard to present a continuum of solutions so we restrict ourselves to the discussion of three representative examples. However, we comment that we are able to obtain similar end states for a wide range of migration parameters. Before presenting these examples we briefly indicate how the librations can be stabilized by the migration process.

As above, the parameters $\langle F_i\rangle $ and $\tau_{\rm c}$ are kept constant, so maintaining D constant, for the duration of the simulation, with $\tau_{\rm c}$ being determined for the initial location of the outer planet. As discussed above, these values may scale with the radius of the planet and we thus expect them to change during migration. However, the semi major axis of the outer planet changes only by $\sim $30% during the simulation consequently we ignore this effect.

 \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig8.eps} \end{figure} Figure 8:

These plots show another possible formation scenario of (Obj)HD 128311(Obj). The damping parameters are the same as in Fig. 7.
Open with DEXTER

5.2 Adiabatic invariance and the stabilization of librations through damping

Both the fast angle $\phi _1$ and the slow angle $\phi _2$ obey stochastically forced oscillator equations. Although we do not have a simple expressions for the oscillation frequencies we know that for fixed eccentricity and planet masses they both scale as the orbital frequency, n1.When the planets migrate inwards together, this increases slowly with time. Accordingly the theory of adiabatic invariance indicates that the libration amplitude should scale as n1-1/2(Peale 1976). If there were no stochastic forcing or other effects, the libration amplitude would decrease on twice the migration time scale $2\tau_a$ once steady eccentricities were formed. This provides only a small effect given the typical amount of migration $\sim $30% in our simulations.

A comparison of the behaviour of a model of the (Obj)GJ 876(Obj) system undergoing convergent orbital migration resulting in the fomation of a commensurabilty both with and without eccentricity damping was undertaken by Lee & Peale (2001). Small libration amplitudes at a similar level are attained in both models. The case without eccentricity damping attains higher eccentricities and somewhat larger libration amplitudes as the evolution continues whereas a low steady libration amplitude is maintained in the case with eccentricity damping. The above situation is consistent with the initial reduction in libration amplitude that occurs prior to the attainment of steady eccentricities as not being due to eccentricity damping but the operation of adiabatic invariance (Peale 1976).

However, a study of simulations of the (Obj)GJ 876(Obj) system without stochastic forcing reveals that eccentricty damping does play a role in damping the slow mode (associated with apsidal corotation) on the eccentricity damping timescale, whereas the fast mode appears to be only marginally affected at least at small amplitudes. Thus the residual librations in these cases consist mainly if not entirely of the fast mode. Therefore the inclusion of eccentricity damping does impart some stability to the commensurability through control of the slow mode and thus its possible interactions with the fast mode. As the slow mode involves $\zeta $ and e1,(see Fig. 2), it is not surprising that this is affected by eccentricity damping.

Note too that because the migration process itself causes libration amplitude reduction on formation of the commensurability some stability is also implied towards some types of disturbance that would take the system away from resonance.

5.3 Model 1

The planets are initially in circular orbits at radii a1=4 AU and a2=2 AU. Stochastic forcing is applied to the outer planet only with the diffusion coefficient $D=6.4\times 10^{-3}~ {\rm cm}^2/{\rm s}^3.$Note that although this value is 640 times larger than the scaling given by Eq. (65), corresponding to a force that is 25 times larger than the simple estimate given in Sect. 3.1, corresponding to smaller reduction factors resulting from a gap or dead zone, it is too small for the angle $\phi _1$ to be brought to circulation during our runs. The results given above and summarized in Fig. 5, and other tests, indicate that similar results would be obtained if D is reduced, but on a longer time scale ${\propto} 1/D,$ provided the migration rate is also appropriately ultimately reduced so that the system can survive for long enough to enable $\zeta $ to be driven into circulation.

The outer planet is made to migrate inwards on a timescale $\tau_{a_1}=8000$ years. The inner planet migrated slowly outwards on a timescaly of $\tau_{a_2}=-20~000$ years. This results in convergent migration. For both planets we use an eccentricity to semi major axis damping ratio of K=8. The resulting evolutionary timescales are significantly larger than those used by e.g. Sándor & Kley (2006) and more easily justified by hydrodynamical simulations.

The time evolution is shown in the left plot of Fig. 6. After resonance capturing all angles are either initially librating or on the border between libration and circulation. The slow mode has a period of $\sim $700 years, whereas the slow mode period is 20 times shorter.

It can be seen from Fig. 6 and other related figures below, that while migration continues libration ampliudes tend to be controlled apart from when e1 becomes either zero or very close to zero. Then, due to stochastic forcing, $\phi _2$ and $\zeta $ start circulating. Subsequently libration is recovered over time intervals for which e1does not attain very small values but eventually additional stochastic forcing together with the repeated attainment of small values for e1causes $\phi _2$ to remain circulating for the remainder of the simulation. After 13 000 years, both the forces due to migration and turbulence are reduced smoothly on a timescale of 2000 years. The result is a stable configuration that resembles the observed system very well.

5.4 Model 2

For this simulation, we lengthened the migration timescales by a factor of 2 to show that the results are generic. We also decreased the value of K to 5.5. This results in larger eccentricities and the final state better resembled our representation of the observed system. However, it should be kept in mind that the eccentricities are not well constrained by the observations. All other parameters are the same as for model 1. The time evolution is plotted in Fig. 7. The resonant angles librate immediately after the capture into resonance as predicted for sufficiently slow migration. However, in the same way as for model 1 described above, stochastic forces make $\phi _2$ circulate soon afterwards. Once the migration forces and stochastic forces are removed between 20 000 and 22 000 years, the system stays in a stable configuration with no apsidal corotation.

5.5 Model 3

All parameters are the same as for model 2. Accordingly model 3 represents another statistical realization of that case. The time evolution is plotted in Fig. 8. The final state is very close to the boundary between libration and circulation of $\zeta.$As discussed above, large uncertainties in the orbital parameters means that the actual system could be in such a state.

Due to stochastic forces, which may result from local turbulence, we are able to generate a broad spectrum of model systems. Some of them undergo a strong scattering during the migration process. However, the surving systems resemble the observed system very well. Although the orbital parameters are not well constrained yet, we showed that values in the right range are naturally produced by a turbulent disk.

6 Discussion

In this paper, we have presented a self consistent analytic model applicable to either a single planet or two planets in a mean motion resonance subject to external stochastic forcing. The stochastic forces could result from MRI driven turbulence within the protoplanetary disk but our treatment is equally applicable to any other source.

We considered the evolution of a stochasticly forced two planet system that is initially deep inside a MMR (i.e. the two independent resonant angles librate with small amplitude). Stochastic forces cause libration amplitudes to increase in the mean with time until all resonant angles are driven into circulation at which point the commensurability is lost. Often a strong scattering occurs soon afterwards for sytems composed of planets in the Jovian mass range.

We isolated a fast libration mode which is associated with oscillations of the semi-major axes and a slow libration mode which is mostly associated with the motion of the angle between the apsidal lines of the two planets. These modes respond differently to stochastic forcing, the slow mode being more readily converted to circulation than the fast mode. This slow mode is sensitive to the attainment of small eccentricities which cause rapid transitions between libration and circulation. The amplitude of the fast mode grows more regularly in the mean, with the square of the libration amplitude in most cases increasing linearly with time and being proportional to the diffusion coefficient D. Of course this discussion is simplified and there are limitations. For example if the total mass of the system is reduced, the disruption time eventually becomes comparable to the libration period. In that case the averaging process that we used in the derivation is no longer valid and the lifetime no longer scales as 1/D.

The analytic model was compared with numerical simulations which incorporated stochastic forces. Those forces, parameterized by the mean square values of each component in cylindrical ccoordinates and the autocorrelation time, were applied in a continuous manner giving results that could be directly compared with the analytic model. The simulations were in broad agreement with analytic predictions and we presented illustrative examples of the disruption process. We performed the simulations for a large range of diffusion parameters, planet masses and initial eccentricities to verify the scaling law for the commensurability disruption time summarized hereafter.

To summarize our results, recalling that the slow angle is driven into circulation before the fast angle, so that the ultimate lifetime is determined by the time taken for the fast angle to achieve circulation, we determine the lifetime, tf,using Eq. (60) setting t=tf, $(\Delta Q)^2 = (\Delta Q)^2_0 = 4,$together with $\gamma_f =p =1,$ so obtaining

$\displaystyle t_f=\frac {a_1^2 \omega_{lf}^2(\Delta Q)^2_0}{36D}\cdot$     (69)

We showed above that this gives good agreement with our numerical results. Using $D=2\langle F_i^2\rangle \tau_{\rm c},$ we can express this result in terms of the relative magnitude of the stochastic forcing in the form

\begin{displaymath}% t_f = 2.4\times 10^{-4}\left(\frac {a_1^2n_1^4}{\langle F_i... ...rt{q_{\rm GJ}}}{n_1\sqrt{q}}\right)^2\frac{q}{q_{\rm GJ}} P_1, \end{displaymath} (70)

where P1 is the orbital period of the outer planet. Here the first quantity in brackets represents the ratio of the square of the central force per unit mass to the mean square stochastic force per unit mass acting on the outer planet. The other quantities in brackets, scaled to the (Obj)GJ 876(Obj) system are expected to be unity, while the last factor $q/q_{\rm GJ}$ is the ratio of the total mass ratio of the system to the same quantity for (Obj)GJ 876(Obj). Here it is assumed that the two planets in the system have comparable masses.

From (70) we see that a non migrating system such as (Obj)GJ 876(Obj) could survive in resonance for $t_f \sim 10^6$ years if the stochastic force amplitude is ${\sim} 10^{-5}$ times the central force. This expression enables scaling to other systems at other disk locations for other stochastic forcing amplitudes. Inference of survival probabilities for particular systems depends on many uncertain aspects, such as the protoplanetary disk model and the strength of MRI turbulence. However, the mass ratio dependence in Eq. (70) indicates that survival is favored for more massive systems. At the present time the number of observed resonant systems is too small for definitive conclusions to be made. However, the fact that several systems exhibiting commensurabilities have been observed indicates that resonances are not always completely disrupted by stochastic forces due to turbulence, but rather may be modified as in our study of (Obj)HD 128311(Obj).

The (Obj)HD 128311(Obj) system is such that the fast mode librates with the slow mode being near the borderline between libration and circulation. We found that such a configuration was readily produced in a scenario in which the commensurability was formed through a temporary period of convergent migration with the addition of stochastic forcing. During a migration phase moderate adiabatic invariance applied to the libration modes together with eccentricity damping leads to increased stabilization and a longer lifetime for the resonance. However, as discussed above, the time evolution of the eccentricity and in particular the attainment of small values plays an important role in causing the slow mode to circulate corresponding to the loss of apsidal corotation. Thus we expect that a large eccentricity damping rate does not necessarily stabilize the apsidal corotation of thesystem. Additional simulations have shown that this is lost more easily for large damping rates ($K\gg 10$). On the other hand it should be noted that this has less significance when one of the orbits becomes nearly circular.

Further observations of extrasolar planetary systems leading to better statistics may lead to an improved situation for assessing the role of stochastic forcing in future.

Acknowledgements

We thank an anonymous referee for useful comments and suggestions. H.R. was supported by an Isaac Newton studentship, the award of a STFC studentship and St John's College Cambridge. Numerical simulations were performed on the Hyades cluster.

References

     
  • Adams, F. C., Laughlin, G., & Bloch, A. M. 2008, ApJ, 683, 1117 In the text  
  • Brouwer, D., & Clemence, G. M. 1961, Methods of celestial mechanics (New York: Academic Press) In the text  
  • Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325 In the text  
  • de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 In the text  
  • Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113 In the text  
  • Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123 In the text  
  • Gammie, C. F. 1996, ApJ, 457, 355 In the text  
  • Hamilton, D. P. 1994, Icarus, 109, 221 In the text  
  • Kasdin, N. J. 1995, Proc. IEEE, 83 In the text  
  • Lecoanet, D., Adams, F. C., & Bloch, A. M. 2008, ArXiv e-prints In the text  
  • Lee, M. H., & Peale, S. J. 2001, BAAS, 33, 1198 In the text  
  • Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 In the text  
  • Nelson, R. P. 2005, A&A, 443, 1067 In the text  
  • Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 In the text  
  • Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805 In the text  
  • Papaloizou, J. C. B. 2003, Celest. Mech. Dyn. Astron., 87, 53 In the text  
  • Papaloizou, J. C. B., & Szuszkiewicz, E. 2005, MNRAS, 363, 153 In the text  
  • Papaloizou, J. C. B., & Terquem, C. 2006, Rep. Progr. Phys., 69, 119 In the text  
  • Peale, S. J. 1976, ARA&A, 14, 215 In the text  
  • Press, W., et al. 1992, Numerical recipes in C (Cambridge: Cambridge University Press) In the text  
  • Reipurth, B., Jewitt, D., & Keil, K. 2007, Protostars and Planets V In the text  
  • Rivera, E. J., Lissauer, J. J., Butler, R. P., et al. 2005, ApJ, 634, 625  
  • Sándor, Z., & Kley, W. 2006, A&A, 451, L31 In the text  
  • Schneider, J. 2009, http://exoplanet.eu In the text  
  • Sinclair, A. T. 1975, MNRAS, 171, 59 In the text  
  • Snellgrove, M. D., Papaloizou, J. C. B., & Nelson, R. P. 2001, A&A, 374, 1092 In the text  
  • Stoer, J., & Bulirsch, R. 2002, Introduction to Numerical Analysis (Springer) In the text  
  • Udry, S., Fischer, D., & Queloz, D. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 685 In the text  
  • Vogt, S. S., Butler, R. P., Marcy, G. W., et al. 2005, ApJ, 632, 638 In the text  
  • Weidenschilling, S. J. 1977, Ap&SS, 51, 153 In the text

All Tables

Table 1:   Parameters of the model systems considered.

All Figures

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330f1a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f1b.eps}\end{figure} Figure 1:

Time evolution of the longitude of periastron $\varpi $, the eccentricity e and the semi major axis a for a single planet. The initial orbital parameters of the planet were taken to be those of (Obj)GJ 876 b(Obj) (see Table 1). Six different realisations starting from the same initial conditions are shown in each panel. The diffusion coefficients are $D = 8.2 \times 10^{-5} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the upper three panels and $D = 8.2 \times 10^{-3} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the lower three respectively. The solid lines correspond to the analytic predictions for the amount of spreading (see text).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig2.eps} \par\end{figure} Figure 2:

Time evolution of the resonant angles and the eccentricities in the system (Obj)GJ 876(Obj) without turbulent forcing. The dominance of the fast mode with period ${\sim } 1.4$ years in the oscillations of $\phi _1$, and the dominance of the slow mode with period ${\sim } 10$ years in $\zeta $ can be clearly seen.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330f3a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f3b.eps}\end{figure} Figure 3:

Time evolution of the eccentricities in the (Obj)GJ 876(Obj) system with turbulent forcing included. The diffusion coefficient is $D = 8.2 \times 10^{-5} ~{{\rm cm}^2}/{{\rm s}^{3}}$ for the two upper plots and $D = 8.2\times10^{-4} ~ {\rm cm}^2/{\rm s}^{3}$ for the two lower plots. The first (uppermost) and the third plots each show a single run. The second and fourth (lowermost) plots show the time averaged eccentricity for four runs. The averaging interval is 1000 years.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330f4a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f4b.eps}\end{figure} Figure 4:

Time evolution of the resonant angles, the period ratio P2/P1 and the eccentricity, e1, in the (Obj)GJ 876(Obj) system with stochastic forcing corresponding to $D = 0.42~ {\rm cm}^2/{\rm s}^{3}$. The vertical lines indicate when the angles enter circulation for a prolonged period. The realisation illustrated in the lower panel scatters shortly after $\phi _1$ goes into circularization.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330f5a.eps}\par\vspace*{1.2mm} \includegraphics[width=9cm,clip]{1330f5b.eps} \end{figure} Figure 5:

Average time until circulation of the resonant angles $\phi _1$ ( top) and $\phi _2$ ( bottom) in the GJ 876 (indicated with +) , GJ 876LM (indicated with $\times $), GJ 876SE (indicated with $\star $) and GJ 876E (indicated with $\Box $) systems as a function of the stochastic forcing diffusion coefficient D.In the case of the upper panel, the analytic curves explained in the text, are from top to bottom for the GJ 876, GJ 876LM, GJ 876SE and GJ 876E systems respectively. In the lower panel the analytic curves, as explained in the text, are from top to bottom for the GJ 876 (solid curve and top dashed curve), GJ 876LM, GJ 876SE and GJ 876E systems respectively.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig6.eps} \end{figure} Figure 6:

The plots on the left hand side show a possible formation scenario of (Obj)HD 128311(Obj). We plot the observed system on the right hand side as a comparison (see Table 1 and text). The plots show the resonant angles $\phi _1, \phi _2, \zeta $, the eccentricities e1, e2 and period ratio P1/P2 for formation scenarios including turbulence and migration. Resonance capturing occurs after 2000 years.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig7.eps} \end{figure} Figure 7:

The plots on the left hand side show another possible formation scenario of (Obj)HD 128311(Obj). Again, we plot the observed system on the right hand side as a comparison (see Table 1 and text). The migration timescales in this run are $\tau_{a,1}=16~000$, $\tau_{a,2}=-40~000$ and K=5.5. Note that the eccentricities are larger compared to Fig. 6 because K is smaller.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{1330fig8.eps} \end{figure} Figure 8:

These plots show another possible formation scenario of (Obj)HD 128311(Obj). The damping parameters are the same as in Fig. 7.
Open with DEXTER
In the text

Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.