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

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
of order
,
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
which is often characterized using the Shakura & Sunyaev (1973)
parameter, with
being the local sound speed and
being
the orbital period. For the most likely situation of small or zero net magnetic flux,
recent MHD simulations have indicated
.
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)

=

(1)

=

(2)

=

(3)

=

(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

(5)

(6)

where
is the gravitational constant, a the semi-major axis and e the eccentricity.
The mean longtitude is
where n the mean motion, with
t0 being the time of periastron passage and
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
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 .
Because of this we may determine them by considering forces of the form
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
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

(7)

The additional terms proportional to the components of
correspond to the Gaussian
form of the equations of motion (Brouwer & Clemence 1961).

The various derivatives involving
can be calculated by elementary means and expressed in terms of
and
One thus finds additional contributions to the equations of motion (1)-(4), indicated with a subscript F, in the form

=

(8)

=

(9)

=

(10)

=

(11)

=

(12)

=

(13)

where the true anomaly f is defined as the difference between the true longitude and the longitude of periastron,
Note that from (9)
we obtain

We remark that the above formalism results
in Eq. (11) which
gives an expression for
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
and
as
dynamical variables rather than e and
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 1/e do not appear.
This can be seen from (11) and (16) which give
in the small e limit

(17)

(18)

Abrupt changes to
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
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
of the inner planet of reduced mass m2is measured from Mand that of the outer planet,
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

(19)

Here
M1=M+m1,
M2= M + m2 and
.
The Hamiltonian can be expressed in terms of
,
i=1,2
and the time t.
The energies are given by
and the angular momenta
with
ai and ei denoting the semi-major axes and eccentricities
respectively. The mean motions are

The Hamiltonian may quite generally
be expanded in a Fourier series
involving linear combinations of the three angular differences
and
(e.g. Brouwer & Clemence 1961).

Near a first order p+1: p resonance, we expect that both
and
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
and 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+ H12, where

(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:

=

(21)

=

(22)

=

(23)

=

(24)

(25)

(26)

Here

(27)

(28)

and

(29)

Note that
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
and
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
and the slow mode has libration frequency
.
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
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
and
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

(30)

where

(31)

and we have used the resonance condition that
(p+1)n1 = p2n2which is satisfied to within a correction of order
Note that for this mode the fact that
very nearly, implies that
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
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
where
Subtracting Eq. (26) from Eq. (25),
differentiating with respect
to time and using this condition results in an equation for

(32)

where

(33)

with

and

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

(34)

The quantities
are readily obtained for each planet from (11).
From this we see that for small eccentricities,
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
in the form

(35)

If equations for quantities that are regular functions of
the Cartesian like coordinates
and
are found, the source terms arising from the external forcing are regular as
.
To illustrate this we consider the pair

It is readily seen that their time derivatives are given by

(36)

(37)

here |U denotes evaluation without external forcing.
There is now no divergence of the external forcing terms as
If as above small amplitude slow mode librations of
about zero are considered,
given that without external forcing
the first term on the right hand side
of (36) is either quadratic in
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
may be replaced by
unity. On taking the time derivatives of (36) and (37),
we obtain under the same approximation scheme

(38)

(39)

Making use
of (32) and its counterparts for
and
for the unforced motion we obtain

(40)

(41)

where
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
or y the slow mode, so we do not
expect such a change of basis to significantly affect conclusions.

We further comment that because
contains
but not
for small eccentricities there are only potential forcing terms 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
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
where the autocorrelation function g(x) is such that
where
is the correlation time and
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
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
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

=

=

(42)

here we take the limit where
is very large corresponding
to an integration time of very many correlation times and
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

(43)

where we recall that
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
without changing the results given below. One readily finds that
Eq. (42) is replaced by

=

=

(44)

where

Note that when
corresponding
to the correlation time being much less than the orbital
period,
For larger
gives a reduction factor
for the amount of stochastic diffusion. For example if we adopt
an exponential form for the autocorrelation function such that

we find

(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

(46)

(47)

(48)

(49)

(50)

We note again that the 1/e2 dependence of
which arises from the coordinate singularity discussed in Sect. 2.1
does not appear for
and
Note that the definition of (h,k) imply consistently with the above that

(51)

(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

(53)

where
There are corresponding expressions that can be obtained from
Eqs. (41) and (35)
for
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

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

=

=

(57)

We then find, working
in the limit of a small eccentricity e1 (but not necessesarily e2), from the corresponding equation for
that

In dealing with Eq. (56) we neglect
in S because after integration by parts
this leads to a contribution on the order
smaller than that derived from
Thus we simply
obtain

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

As we are interested in the case p=1, the width of the libration zone
is
we see that the time for
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
or
as that
obtained for
or
from Eq. (49)
applied to the outer planet.
This indicates that ,
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
for small e1.
The latter quantity was assumed constant in the analysis.
As implied by the discussion of Sect. 2.1
abrupt changes to
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
,
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
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
is significantly smaller than Q unless e1 starts out
with a very small value. This would suggest that Q reaches circulation
before
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
remains librating as it was initially.
One expects to recover standard
secular dynamics for
from the governing Eqs. (21)-(25),
when these are averaged over an assumed
circulating
with constant
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
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 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
which induce
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)
and the auto correlation time
.

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

For planets under the gravitational influence
of a protoplanetary disk, the natural scale for the force per unit mass components, Fi, is
,
where
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
at a distance
from its centre assuming
that all its mass is concentrated there. The result is independent of
The natural correlation time
is the inverse of the orbital angular frequency
.
To set the natural scale for D, we adopt a
minimum mass solar nebula model (MMSN, see Weidenschilling 1977) with

(63)

This provides a natural scale for D as a function of the local disk radius and the central stellar mass
through

(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
(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
seems reasonable although it might be even smaller (de Val-Borro et al. 2006).
The correlation time
is actually found to be
approximately
(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
and we expect a natural scale for the diffusion coefficient to be
specified through

(65)

The same value of D may be equivalently scaled to the orbital parameters
of the planets without reference to the disk by writing

(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
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
(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
.
In this simulation, we use the reduced diffusion coefficients

D

=

(67)

D

=

(68)

which can be represented by a correlation time of half the orbital period
and a specific force with root mean square values
and
,
respectively
(see also Eq. (66) above).

The resulting random walks undergone by e, a and
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.

Figure 1:

Time evolution of the longitude of periastron ,
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
for the upper three panels
and
for
the lower three respectively.
The solid lines correspond to the analytic predictions for the amount of spreading
(see text).

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
while also being present in those of
On the other hand the slow mode, which is seen to have a
period of about 10 years,
dominates the librations of
while also being present
in those of e1 and

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

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
years in the oscillations
of ,
and the dominance of the slow mode
with period years in
can be clearly seen.

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 then the period of the fast mode scales as
and the period of the slow mode scales as

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
and
.
We remark that all of our simulations were run with constant
values of
This was done maintaining the parameters
and
to be constant, with
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.

Figure 3:

Time evolution of the eccentricities
in the (Obj)GJ 876(Obj) system with turbulent forcing
included. The diffusion coefficient is
for the two upper plots and
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.

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
to undergo large oscillations
and eventually circulation (see Eq. (49)).
Should stochastic forces cause
e1 to reach zero, on account of the coordinate singularity
becomes undefined.
Subsequently very small perturbations are able to produce a
small eccentricity
with
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
with values
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
enters circulation
prior to the fast angle
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
enters circulation, stochastic forcing
acts to increase the libration amplitude of the fast mode.
This mode dominates both
the librations of the resonant angle
and the semi major axes.
Eventually
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
In this context we note that reducing D increases
the evolutionary time which has been found, both analyticaly and numerically
to be
(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.

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
.
The vertical lines indicate when the angles enter circulation
for a prolonged period.
The realisation illustrated in the lower panel
scatters shortly after
goes
into circularization.

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
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
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
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 We also see from Fig. 4
that the angle
which has a large contribution from the slow mode,
behaves in the same way as 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 attains circulation.
Also as expected, the time when
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
is always taken
to be given by
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 2 in all our simulations.

In this context we consider the fast angle
and the slow angle though, as we saw above, the latter can be replaced by which exhibits the same behaviour.
As
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
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.

Figure 5:

Average time until circulation of the resonant angles
( top) and
( bottom) in the GJ 876
(indicated with +) , GJ 876LM (indicated with ), GJ 876SE (indicated with )
and GJ 876E (indicated with )
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.

From Fig. 5 it is apparent that the evolutionary times scale is
for D varying by many orders of magnitude. The analytic estimates for the libration
survival times of 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
based on Eq. (59) using the initial value
of e1 overestimates the lifetime by a factor of at least
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
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 .
We estimate the lifetime
of librations of
by calculating the time
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 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
These simple estimates are in good agreement with the numerical results and
accordingly support the idea that
and
are indeed coupled in the non linear regime.
For low mass planets this simple analytic prediction underestimates the lifetime of
by a factor of 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
is unchanged, as this does not depend significantly on the eccentricity.
However, the lifetime of the librating state of
is a factor of 2-3 longer.
This is due to the fact that a larger excitation of
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
librating and the angle
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
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
.

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
,
the eccentricities e1, e2 and period ratio P1/P2 for formation scenarios including turbulence and migration.
Resonance capturing occurs after 2000 years.

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
,
and K=5.5. Note that the eccentricities are larger compared to Fig. 6 because K is smaller.

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,
and
.
We also define the ratio of the above timescales to be
.
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
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
and
are kept constant, so maintaining
D constant, for the duration of the simulation, with
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 30% during the simulation
consequently we ignore this effect.

Figure 8:

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

5.2 Adiabatic invariance and the stabilization of librations through damping

Both the fast angle
and the slow angle
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
once steady eccentricities were formed.
This provides only a small effect given the typical amount
of migration 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
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
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
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
provided the migration rate is also appropriately ultimately reduced
so that the system can survive for long enough to enable
to be driven
into circulation.

The outer planet is made to migrate inwards on a timescale
years. The inner planet migrated slowly outwards on a timescaly of
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 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,
and
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
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
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 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,
together with
so obtaining

(69)

We showed above that this gives good agreement with our numerical results.
Using
we can express this result in terms
of the relative magnitude of the stochastic forcing in the form

(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
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
years if the stochastic force amplitude
is
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 ().
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.

All Tables

All Figures

Figure 1:

Time evolution of the longitude of periastron ,
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
for the upper three panels
and
for
the lower three respectively.
The solid lines correspond to the analytic predictions for the amount of spreading
(see text).

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
years in the oscillations
of ,
and the dominance of the slow mode
with period years in
can be clearly seen.

Time evolution of the eccentricities
in the (Obj)GJ 876(Obj) system with turbulent forcing
included. The diffusion coefficient is
for the two upper plots and
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.

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
.
The vertical lines indicate when the angles enter circulation
for a prolonged period.
The realisation illustrated in the lower panel
scatters shortly after
goes
into circularization.

Average time until circulation of the resonant angles
( top) and
( bottom) in the GJ 876
(indicated with +) , GJ 876LM (indicated with ), GJ 876SE (indicated with )
and GJ 876E (indicated with )
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.

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
,
the eccentricities e1, e2 and period ratio P1/P2 for formation scenarios including turbulence and migration.
Resonance capturing occurs after 2000 years.

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
,
and K=5.5. Note that the eccentricities are larger compared to Fig. 6 because K is smaller.

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.