A dynamical analysis of the Kepler-11 planetary system

# A dynamical analysis of the Kepler-11 planetary system

Cezary Migaszewski, Mariusz Słonina and Krzysztof Goździewski11footnotemark: 1
Toruń Centre for Astronomy, Nicolaus Copernicus University, Gagarin Str. 11, 87-100 Toruń, Poland
E-mail: c.migaszewski@astri.umk.plE-mail: m.slonina@astri.umk.plE-mail: k.gozdziewski@astri.umk.pl
Accepted 2012 August 23. Received 2012 August 20; in original form 2012 May 3
###### Abstract

The Kepler-11 star hosts at least six transiting super-Earth planets detected through the precise photometric observations of the Kepler mission (Lissauer et al.). In this paper, we re-analyze the available Kepler data, using the direct -body approach rather than an indirect TTV method in the discovery paper. The orbital modeling in the realm of the direct approach relies on the whole data set, not only on the mid–transits times. Most of the results in the original paper are confirmed and extended. We constrained the mass of the outermost planet g to less than 30 Earth masses. The mutual inclinations between orbits b and c as well as between orbits d and e are determined with a good precision, in the range of [1,5] degrees. Having several solutions to four qualitative orbital models of the Kepler-11 system, we analyze its global dynamics with the help of dynamical maps. They reveal a sophisticated structure of the phase space, with narrow regions of regular motion. The dynamics are governed by a dense net of three– and four–body mean motion resonances, forming the Arnold web. Overlapping of these resonances is a main source of instability. We found that the Kepler-11 system may be long-term stable only in particular multiple resonant configurations with small relative inclinations. The mass-radius data derived for all companions reveal a clear anti-correlation between the mean density of the planets with their distance from the star. This may reflect the formation and early evolution history of the system.

###### keywords:
celestial mechanics – planetary transits – Kepler-11 – Arnold web
pagerange: A dynamical analysis of the Kepler-11 planetary systemReferencespubyear: 2012

## 1 Introduction

The Kepler space mission is a breakthrough in the field of searches for the Earth–like extrasolar planets (Borucki et al., 2010; Koch et al., 2010; Jenkins et al., 2010; Caldwell et al., 2010). About of 150,000 solar dwarfs are monitored by 0.95—meter Kepler telescope. The photometric data are publicly available from the MAST archive.

To date, the mission identified more that 2,200 planetary candidates (Batalha et al., 2012). Among them, many multi-planet systems are found. For instance, planets were confirmed in two–planet configurations, i.e., Kepler-10 (Batalha et al., 2011; Fressin et al., 2011), Kepler-25, 26, 27, 28 (Steffen et al., 2012), Kepler-29, 31, 32 (Fabrycky et al., 2012), Kepler-23, 24 (Ford et al., 2012); in three-planet systems Kepler-9 (Holman et al., 2010), Kepler-30 (Fabrycky et al., 2012), Kepler-18 (Cochran et al., 2011); in four–planet systems (Borucki et al., 2011), as well as in five–planet configurations Kepler-20 (Gautier et al., 2011; Fressin et al., 2011), Kepler-33 (Lissauer et al., 2012). The Kepler-11 hosts six planetary companions (Lissauer et al., 2011). The transiting planet candidates can be confirmed through determining their masses with the help of the so-called Transit Timing Variations method (TTV, Holman & Murray, 2005; Agol et al., 2005). In this approach, the (O-C) variations between observed mid-transit times and their ephemeris are the observables, which can be fitted by an appropriate orbital model. In recent papers, also additional observables are analysed, like the so-called Transits Duration Variations (TDVs) (see, e.g., Nesvorný et al., 2012).

In this paper, we re-analyse the photometric data of Kepler-11 with a modified, direct approach providing an alternate estimation of masses and orbital elements. To describe this method further in the paper, we recall shortly the main conclusions in (Lissauer et al., 2011). Using the TTV method and an assumption of strictly coplanar model of the system, they determined masses of five inner planets in the range of a few Earth masses. The outermost planet interacts weakly with the inner companions, and its mass could be roughly constrained as smaller than the Jupiter mass. It has been not confirmed as a planet, although the probability of blending is very small, . Orbital eccentricities in the Kepler-11 system were determined only for the five inner objects. Due to the assumption of coplanarity, a determination of mutual inclinations between the orbits was not possible. Lissauer et al. (2011) argue that these inclinations should remain in the range of [0,2] degrees. The dynamical analysis have revealed that the system is not involved in the mean motion resonances (MMRs), however a pair of planet b and planet c is close to 5:4 MMR.

The determined masses and radii of the planets imply constrains on their chemical composition. Planets d, e and f might have similar internal compositions to those of Uranus or Neptune, while planets b and c are rather ice rich, with a smaller amount of H/He mixture than these planets in the Solar system.

In this paper, we focus mostly on the global dynamics of the system and a few aspects which were not addressed in the discovery paper.

First of all, we model the available Kepler data through a direct algorithm that relies on the self-consistent -body fitting of the light-curves, instead of the TTV method applied in the discovery work. The TTV algorithm makes use of the transit times a posteriori, after they are determined from the light curves. Through extensive numerical experiments, we found that the direct approach brings more information than the TTV method. For instance, we could constrain the mass of the outermost planet to less than  Earth masses. We also found significant bounds for the mutual inclinations to less than for planets b and c as well as for planets d and e.

The direct model, also called the dynamical-photometric model, already was used in a few papers. For instance, it was applied to analyse the light curve of the triple-star system KOI-126 (Carter et al., 2011), and to estimate masses of two planets transiting Kepler-36 (Carter et al., 2012). This algorithm also verified the Kepler-9 model, which was found first with the help of the TTV algorithm (Holman et al., 2010).

A number of initial conditions found with the direct approach makes it possible to investigate the dynamics of the system. We focus on the short–term time scale, governed by the mean motion resonances. We study the multi-dimensional structure of the phase space with the help of dynamical maps. In the vicinity a few qualitative transit models considered in this work, the dynamics are governed by a dense net of 3–body and 4–body mean motion resonances. This net may be identified with the Arnold web, which is a feature of close to integrable Hamiltonian systems. The Kepler-11 appears as strongly resonant extrasolar system, and this feature may reflect its trapping into MMRs at the early stages of the formation and evolution.

Using a new determination of the masses and radii, we found a curious mass-radius relation implying a clear anti-correlation between the mean density of the planets and their distances from the star. Their densities exhibit a sequence of planet b which is denser than Neptune, through the Neptune-like planet c, Uranus-like planet d, Jupiter-like planets e and f, and planet g which is likely Jupiter/Saturn-like.

The paper is structured as follows. In Sect. 2, we shortly describe the photometric data of Kepler-11 available in the MAST archive. We also refine the observational TTV model. In Sect. 3, we present the results derived through intensive computations with the bootstrap algorithm. Furthermore, we discuss a possible composition of the planets (Sect. 4). Section 5 is devoted to the dynamical analysis of the Kepler-11 system. Conclusions and prospects for a future work are given at the end of this paper.

## 2 Transits in a multi-planet system

The photometric data of Kepler-11 were taken from the MAST archive. At the time of writing this paper, the publicly available light-curves span about of  days in six parts. These data were binned on -minute intervals. We analysed a “de-trended” data set derived through a smoothing procedure. At first, we isolate all transits from the light curve. Then the moving average with a time-step of  days provides the mean level of the flux. Next, we construct an interpolated, reference light curve with the cubic spline on these nodes. Finally, we divide the raw flux, with all transits data, by its values of the reference, mean level flux curve.

The de-trended data available in the MAST database exhibit a growth of the flux shortly before and after a particular transit. In some parts of the available light-curves, spanning approximately  days, the measurements appear in the raw form. We did not use these data, aiming to analyze a possibly uniform set of observations.

### 2.1 Modeling the stellar flux

A common model of photometric observations of a star transited by planetary companions consists of two major parts. The first part concerns the flux deficit due to small, dark objects passing in front of the star. At first, the average orbital periods are determined. Then transit depths and duration times are parametrized on the basis of phase–folded light-curves. Single mid-transit times are also determined. At the second level, we can estimate the planetary masses and orbital elements fitting a model of motion of mutually interacting planets.

We focus on the first level of the photometric analysis. To compute the flux deficit, we use the quadratic limb darkening model (Mandel & Agol, 2002), recalling that the Kepler-11 light-curves are relatively noisy and sampled with a low frequency,

 ΔI(r)=1−γ1(1−cosθ)−γ2(1−cosθ)2, (1)

where is the normalized radial coordinate w.r.t. the centre of the stellar disk, is the angle between the direction to the observer and the normal to the stellar surface. The two limb-darkening coefficients and must be positive and (see a study of the limb darkening coefficients for a few target stars of the Kepler mission, Howarth, 2011). For small ratio of planet radius to the stellar radius , Mandel & Agol (2002) found an analytic approximation of the flux deficit, which depends on the normalized distance  between the centers of stellar and planetary disks, projected onto the sky plane (see Eq. 8 in the cited paper), as well as on and .

If more than one planet transits the star at the same time, the total flux deficit can be computed as the sum of the deficits caused by particular planets. Obviously, are the same for all planets, while and are different for each object. If transiting planets are small, we can use a simple model of independent transits rather than more general treatment (e.g., Pál, 2011). Because we model the photometric measurements directly, by reconstructing the whole light-curve, we are not restricted to single transits and mid-transit times. Also multiple transits can be covered. In light of relatively narrow observational window, multiple transits are very helpful to constrain orbital elements of the transit model.

Figure 1 displays a few selected fragments of the data set marked with red dots and error bars which are over-plotted on the synthetic curve best–fitting the data (blue curve). The fitting procedure will be described in more detail in Sect. 2.3. The last panel shows transits of three planets (b, d and e).

### 2.2 The model of orbital motion

The orbital motion of multiple planetary system is described in terms of the full -body problem in the Poincaré reference frame (e.g., Morbidelli, 2002). In this frame, the Cartesian coordinates of the planets are astrocentric, while their velocities are barycentric. The equations of motion are integrated with the second order symplectic integrator SABA (Wisdom & Holman, 1991; Laskar & Robutel, 2001). It provides 2-3 times better CPU performance than other algorithms, which we tested (like the Bulirsh-Stoer-Gragg scheme, BGS) constrained with the same time–step accuracy. To speed-up the computations even more, we did not integrate the system at all measurements moments. This would force -minute step-size of the integrator. Instead, we fixed this step-size to of the innermost period of planet b, i.e.,  day. Furthermore, the flux function is computed only close to the mid-transits. Ingress and egress times of particular events are tabulated. When a transit takes place, the coordinates of particular planet at time required to evaluate the flux deficit are determined through the polynomial interpolation on five nodes around . Through a comparison with the direct, full-accuracy integrations with the BGS algorithm, we found that the selected time time-step and the number of interpolation nodes provide a sufficient precision and acceptable CPU overhead. We examined this method by changing the number of nodes in the polynomial interpolation, as well as the time step-size. The flux level, interpolated on five nodes and with  day, differs from its exact value by less than .

### 2.3 Optimization algorithm and error estimation

We searched for the best–fit model of the transits by a common minimization of the function. This function is defined as follows:

 χ2ν=1N\scriptsize obs−N\scriptsize p−1N\scriptsize obs∑j=11σ2j[Fj−F(tj)]2, (2)

where is the number of observations, is the number of free parameters, is the number of the degrees of freedom, is the error of the -th observation , and is a model function evaluated at time . This form of the –function is correct if the uncertainties are uncorrelated (see, e.g., Baluev, 2009). To verify whether the available photometric data fulfill this assumption, one has to use a more general statistical model incorporating the red–noise effect. However, under particular settings of our -body photometric model, this would require an enormous CPU overhead. Hence, we use equation 2 as a reasonable first order approximation.

The best–fit parameters of the transits model are searched through a two–step optimization method. In the first step, we apply a robust and well tested quasi-global Genetic Algorithm (GA, see, Charbonneau, 1995; Deb, 2004)222We use publicly available implementation by Kalyanmoy Deb, see http://www.iitk.ac.in/kangal/pub.htm which makes it possible to find promising solutions. A local, fast gradient method (here, the Levenberg-Marquardt algorithm) is then used to refine the solutions found in the GA step. Such an approach is called the hybrid optimization (see Goździewski et al., 2008, and references therein). Let us note that the parameter space is huge as it has dimension of 50. Some of these parameters can be determined very well, like the orbital periods of the transiting planets. Unfortunately, due to the relatively short observational time window, many parameters which are critical for the stability (relative inclinations, masses, nodal lines) cannot be well constrained. It makes the fitting process a challenging problem.

The parameter errors are estimated through the bootstrap algorithm (see, e.g., Press et al., 1992). The bootstrap is CPU–demanding, but it is straightforward method to estimate standard errors in high-dimensional problems and for large number of data. The light–curves which we analyzed have 22,000 points. The bootstrap algorithm requires to find the best–fit solutions to a large number of synthetic sets derived through random sampling with replacement from the original measurements. To obtain reliable error estimates of the best–fit parameters, one needs at least synthetic solutions. When such a large set of the best–fit models is gathered, we constructed normalized histograms for each free parameter. These histograms reflect the parameter distribution in response to the errors of the measurements, and may be smoothly approximated by an asymmetric Gaussian function. This makes it possible to determine the standard uncertainties. To perform the bootstrap procedure, at first one needs to find reliable best–fit parameters for the nominal data set. This step was done through an intensive quasi-global search with the help of the hybrid algorithm. The bootstrap computations are CPU-time consuming and were performed on the reef CPU-cluster of the Poznań Supercomputing Centre.

### 2.4 Numerical setup of the dynamical analysis

In spite of small eccentricities and apparently co-planar orbits, the Kepler-11 systems is orbitally very active. It appears as dynamically packed planetary system (the definition is given in Barnes et al., 2008), with only narrow stable zones in the phase space. For this reason, we used the best-fit model solutions gathered in the bootstrap search as the input data to extensive dynamical study of this system. As we will discuss later, a study of the stability is a challenging problem. Due to relatively short observational window ( days), weak transits having depths comparable with the measurements errors and a small number of data points covering particular transits (typically ), the derived initial conditions may be shifted away from the real configurations.

To investigate the dynamics of the Kepler-11 system in a global manner, we applied an approach in our previous papers which is well established in the literature. It relies on reconstructing the structure of the phase space with the fast indicator MEGNO (Cincotta & Simó, 2000; Cincotta et al., 2003). This dynamical characteristic makes it possible to distinguish between regular (stable) and irregular (chaotic, unstable) trajectories in the phase space by computing relatively short numerical orbits. Having representative solutions selected in the bootstrap statistics, we study their neighborhood on the dynamical maps. Constructing a dynamical map relies on two model parameters, e.g., the semi-major axes of a pair of planets. The selected parameters are varied in the given range at a discrete grid. The remaining components of the initial parameter vector are fixed at their nominal values. If it is necessary, they are altered to preserve the observational constraints. Then we calculate MEGNO at each point of the grid. Dynamical maps are informative and become a standard numerical tool helpful to understand the global dynamics of multiple systems.

To compute the MEGNO indicator, we must solve the variational equations to the equations of motion of the planetary -body problem. The Kepler-11 system architecture with low–eccentric orbit and small masses is an ideal target for an efficient symplectic algorithm described in (Goździewski, 2003; Goździewski et al., 2008). The general-purpose integrators, like the Runge-Kutta or Bulirsh-Stoer-Gragg schemes are not efficient nor accurate enough in this case. These methods introduce a systematic drift of the energy and other integrals. To avoid such errors, and to solve the variational equations, we apply the tangent map introduced by Mikkola & Innanen (1999). As the very basic step, it requires to differentiate the “drift” and “kick” maps of the standard leap–frog algorithm. The variations may be then propagated within the same symplectic scheme, as the equations of motion. Having the variational vector computed at discrete times, we find temporal and mean of the MEGNO at the -th integrator step , (Cincotta et al., 2003; Goździewski et al., 2008):

 Y(j) = (j−1)Y(j−1)+y(j)j, y(j) = j−1jy(j−1)+2ln(δjδj−1)

with initial conditions , , . The MEGNO maps tend asymptotically to

 Y(j)=ahj+b,

where for quasi-periodic orbits, for stable, periodic orbit, and for chaotic orbit with the maximal Lyapunov exponent . The tangent MEGNO map is linear, hence the variational vector can be normalized, if its value grows too large for chaotic orbits. In practice, we stop the integration if the MEGNO indicator reaches a given limit (usually, ).

The symplectic maps were propagated with the 4-th order SABA scheme in (Laskar & Robutel, 2001). A choice of the fixed step-size must be carefully controlled. We did this, checking whether the relative energy error is “flat” across the dynamical map (Goździewski et al., 2008) and sufficiently small. Indeed, the step–size  day preserved this error at a level of over the total integration times up to  yr ( periods of the outermost planet). This time scale is long enough to detect the most significant 2-body and 3–body MMRs though even such integration period may be insufficient to detect all “dangerous” unstable resonances. Weakly chaotic motions due to multi-body MMRs still may lead to catastrophic events after much longer time (Goździewski et al., 2008).

The dynamical maps in this paper have typical resolution up to pixels. This requires an enormous CPU-time. It is basically not possible to perform such intensive computations on a single workstation. Therefore, we used our new Message Passing Interface (MPI) based environment Mechanic (Slonina et al., 2012) to perform the computations in a reasonable time in CPU-clusters333Informations on this project may found at http://git.astri.umk.pl/projects/mechanic.. They were performed on the reef cluster at the Poznań Supercomputing Centre. A Mechanic run of a typical dynamical map occupied up to 1200 CPU cores for  hours.

### 2.5 Free parameters of the transit model

The free parameters of the transit model are the stellar radius , the limb darkening coefficients ; the mass , radius and orbital elements of each planet in the system, where b,c,d,e,f,g. Planetary orbits are described through the Poincaré geometric, osculating elements at the epoch of the first observation JD 2455964.51128: a tuple is for the semi-major axis, eccentricity, inclination to the plane of the sky, the longitude of ascending nodes, the argument of pericenter, and the mean anomaly, respectively. The orbital node of the first planet, due to invariance of the model with respect to a rotation of the whole system. The inclinations are obviously close to . A deviation from is irrelevant only for single-planet systems.

In a multi-planet system, some orbits may be inclined to the sky plane by angle , which implies different relative inclinations between orbits of particular planets, even for the same longitudes of nodes. Due to the invariance of transits with respect to the direction of the total angular momentum of the system, a combination of means the same geometry as . Thus, when necessary, for a given planet p we can fix the range of and are corrected for remaining companions, in accord with the invariance relation.

Orbital elements are not fully suitable for transiting systems with small relative inclinations and small eccentricities. To avoid singularities and weakly constrained elements, like when (circular or weakly eccentric orbits), we use the Poincaré modified elements , instead of ).

Similarly, the orbital period is more suitable for model fitting than since the semi–major axis depends on the planetary mass (a free parameter) and on the stellar mass , which is fixed to , but it can be also fitted. Hence, we define as one of osculating elements related to through the IIIrd Keplerian law.

The mean anomaly strongly depends on . It determines the relative orbital phase (the mean longitude) but is also related to . This may be avoided by choosing the time of the first transit as a free parameter instead of , because it is one of the directly determined observables from the light-curves. Simple relations between and Poincaré canonical elements may be derived easily.

### 2.6 Direct and indirect transit model parameters

The direct parameters of the transit model are determined from the basic observables: it is the mean period of transits , the depths, duration times of the transits, and shapes of the light-curves (through the limb–darkening coefficients). These data are usually derived from the period–phased light–curves of particular planets. The depths and durations of transit determine the ratio of planetary and stellar radii, . If the stellar mass is fixed then and may be resolved. We can also determine up to the angular momentum direction invariance, and . In general, the mean period of transit events is different from the osculating orbital period at the epoch of the first observation, . A shape of the event–period phased light–curves make it possible to fit the limb darkening coefficients, .

These parameters of the transit model are independent on the the -planet dynamics. Hence the remaining are indirect parameters. To resolve them, a dynamical model of the orbital evolution is required. The indirect parameters consist of planetary masses as well as orbital elements, , , and (instead of ). Knowing and , we may fix the osculating semi-major axis at the date of the first (or prescribed) observation.

We would like to note, that the above distinction for two types of model parameters is somehow arbitrary in our photometric model. In our algorithm both the direct and indirect parameters are fitted simultaneously, unlike, for instance, the TTV algorithm, in which the direct parameters are fitted at the first stage, and the indirect parameters are fitted in the next step.

Usually, the direct parameters can be estimated much more reliably than the indirect parameters. Even a potential derivation of the indirect parameters depend on the particular model of motion. i.e., kinematic — Keplerian, or dynamic — Newtonian, and on the used method of modelling the observations. In the Keplerian (kinematic) model (see, e.g., Agol et al., 2005), mid-transits of a given planet are governed by geometric reflex motion of the star around the center of mass in a sub-system composed of the star and all inner planets. For instance, transits of planet d are affected by planets b and c, but any outer planet does not affect transits of its inner companions. Hence, in accord with the Keplerian model, the indirect parameters of the outermost planet g in the Kepler-11 system cannot be determined at all.

In a given pair of planets, the outer companion affects the transits times of the inner planet only through gravitational mutual perturbations which lead to changes of osculating orbital elements. To account for the mutual interactions, one has to apply the self-consistent -body model of motion of the system.

Usually, to resolve the indirect parameters from photometric observations, the well known TTV method is used (Agol et al., 2005). It has two steps. At first, we determine the mean periods, the mid-transits, and then the (O-C) residua, i.e., differences between the measured and ephemeris transit times. The (O-C) variations are observables in the second step during which we search for masses, eccentricities, and arguments of pericenters of planetary companions. The TTV method in this form has a limitation, because it does not make any use of transit depths nor their duration times. If the individual inclinations of planets are different, the planets transit the parent star usually at different attitudes. Hence the transit depths as well as duration times may vary, like the (O-C) of mid-transits. This information can be used to better constrain the transit model.

The mutual inclinations depend on the longitudes of ascending nodes in accord with

 cosΔIi,j=cosIicosIj−sinIisinIjcos(Ωi−Ωj).

Because the inclinations of transiting planets, (, ) must be close to then . Within this approximation, the TTV method is apparently not sensitive for individual . In fact, different values of imply different mutual inclinations affecting the dynamics and (O-C). However, the dynamical variability of (O-C) due to mutual interactions is weaker than the geometric variability due to changes of transit depths and duration times reflecting the motion of the star around the mass center of the system.

Overall, by direct modeling of the light-curves (photometric measurements), rather than the mid-transit times, we can resolve the (O-C) with an improved precision. Modeling the light-curves in terms of the -body model is CPU-demanding, but it makes it possible to estimate individual longitudes of nodes and mutual inclinations. In particular, as will be shown later, the direct method helped us to derive accurate relative inclinations between planets b and c, as well as between d and e .

## 3 Results of the bootstrap analysis

We performed the direct bootstrap TTV analysis of a few different orbital models of the Kepler-11 system. In the most general case (I), all parameters discussed in the previous section are the free parameters of the fit model. Some of them are poorly constrained by the observations, in particular, the eccentricity of planet g and particular longitudes of nodes. Therefore, we also studied less general models, in which some of weakly constrained parameters are fixed. In the second model (II), , , i.e., . In the third model (III), also , while in the last model (IV), are all fixed at . Because inclinations are not exactly equal to , also .

For each of these four transit models, we applied the bootstrap algorithm and we gathered sets of solutions for each instance of the transit model.

### 3.1 Model I: systems with eg≠0

Figure 2 shows an outcome of the bootstrap algorithm in the form of normalized histograms constructed for and , and depicted from the left to the right panel, respectively. The red solid curves illustrate the best fit asymmetric Gauss function to the histogram bins. The formal errors are marked with red bars displayed above the histograms. The best fit parameters corresponding to the maximum of the Gaussian distribution are written in the respective panels, and they may be compared with the nominal solutions given in Table  1. The uncertainties of the eccentricity and longitude of node of planet g are relatively large.

Because the nominal system is dynamically unstable, we examined the whole set of bootstrap solutions by calculating their MEGNO indicator on the time interval of  yr. It corresponds to periods of the most distant companion. Such a characteristic time scale should be long enough to detect unstable solutions due to low–order 2–body and 3–body mean motion resonances (Goździewski et al., 2008, and references therein). Unfortunately, all initial configurations exhibit large values of , indicating that the system is strongly chaotic. The main source of instability are crossing orbits in the system, that lead to disruptive events , i.e., one or more of the planets were ejected from the system or collided with the parent star. None of the tested solutions passed the direct integration over 10 Myr.

The parameter space of the Kepler-11 system is -dimensional, and the dimension of the phase space of the -body model is -dimensional. The experiments indicate that this system can be locally chaotic and its phase space is filled with mostly unstable solutions. Then only small regions of stable MMRs may be present. In the light of a large dimension of the phase space, the gathered statistics of best–fit configurations is still very poor. We conclude that due to short data span of only  days, and unconstrained elements of the most general model, we cannot find reliably stable solutions assuming the most general transit model I. Unfortunately, in this high–dimensional problem an alternate GAMP algorithm that relies on the optimization with imposed stability constrains (Goździewski et al., 2008) would be CPU-time expensive.

### 3.2 Transit model II: systems with eg=0

In the next model, we narrow the mostly unconstrained parameters of the transit model. We fix the eccentricity , hence and are both equal to . The results of the bootstrap algorithm are illustrated in Figs. 3-9. All panels in these figures are constructed in the same manner as Fig. 2. We tested, whether the best–fit parameters encompass at least marginally stable solutions with after  yr.

Figure 3 shows the normalized histograms for masses of particular planets expressed in the Earth masses. Besides formal uncertainties obtained through the bootstrap (filled red circles), the best-fit parameters derived in (Lissauer et al., 2011) are plotted (blue filled circles). Clearly, these estimates coincide very well in both cases. There is one exception though, since the mass of planet g is not resolved in Lissauer et al. (2011). The direct code helps to resolve also this mass. It is constrained surprisingly well, in spite of a narrow observational window. This result confirms our predictions. Because the orbital model is constrained by all measurements, not the TTVs only, the direct algorithm makes use of dynamical information contained in the transit depths and widths.

For a reference, black and green asterisks in Fig. 3 mark masses of the Uranus and Neptune, respectively. The masses of planets b and f appear in a range specific for the super-Earths. They are significantly smaller than the masses of two most distant planets in our Solar system but, as we will show in the next section, their chemical composition has likely much common with the ice giants in the Solar system.

The next Fig. 4 shows histograms constructed for planetary radii expressed in the unit of the Earth radius. These results confirm data in the discovery paper. Similarly to the previous plots, the radii of Uranus and Neptune are marked with asterisks. They are also labeled with and , respectively. The derived radius of planet g confirms a hypothesis that it may belong to the Uranus/Neptune–class. We note that most of the planets has radii smaller than , and only planet e has its radius larger.

Histograms of the mean densities are presented in Fig.  5. The -axis is for the density expressed w.r.t. the Earth density. Black and green asterisks mark the values characteristic for Uranus and Neptune, respectively. The mean densities of Saturn and water are also marked with the red and blue symbols, respectively. According to this plot, the less dense planet e has a density of Saturn. The most dense planet b may be almost as dense as the Earth. The densities of the other planets span a range characteristic for Saturn and Neptune, from to .

Figure 6 is for the bootstrap histograms constructed of the semi-major axes. These parameters are the best determined among all of the transit models, with uncertainties of the order of only. We do not compare these results with data in (Lissauer et al., 2011) because they accounted for the formal error of the stellar mass. Note that we fixed , because we found that this parameter is unconstrained by the photometric data. Yet it seems that the function monotonically increases in the range of .

The first five panels of Fig. 7 are for the eccentricities, and the bottom, right-hand panel is for . These histograms confirm that the eccentricities of planets b to f are small, typically less that , and the arguments of pericenters are not well constrained. The last panel assures us that is determined with an error of only , recalling a narrow time–window of the photometric data. The best–fit parameters of model II are given in Tab. 2.

Inclination was constrained to the range, and due to the invariance rule implied by the direction of the total angular momentum, the remaining inclination may be smaller and larger than . We tested whether there is a correlation of the transit events with a given half–disc of the star. We found that both cases are equally possible. Because the orbits are inclined to the plane of the sky at angles close to , the relative inclinations with the same longitudes of nodes may be . As expected, the indirect parameters are unconstrained, see Tab. 2. Therefore, the main contribution to the uncertainties of the relative inclinations comes from ambiguous estimates of rather than of .

Curiously, there appears a clear correlation between mutual inclinations in particular pairs of orbits, namely c and e, f and e, as well as d and f. This can be seen in normalized histograms constructed for the inclinations, Fig. 8. For a chosen planet, we transform to range (in accord with the inclination invariance rule), and we compute the bootstrap histogram for . Panels of Fig. 8, from the left to the right, are for pairs . If then much more likely than . Similarly, if , then appears more likely than .

For particular pairs of planets, the relative inclinations can be determined surprisingly well. Figure 9 shows the bootstrap histograms for such pairs which exhibit well constrained values. These histograms reveal that orbits of planets b and c are almost coplanar. Similarly, the pair of planets d and e form an almost coplanar sub-system. The mutual inclinations of orbits in these pairs are less than , with most likely values of . The remaining panels indicate that the mutual inclinations between five inner orbits remain within a few degrees range. Their upper limits are not so small as in the first two sub-systems. The outermost orbit of planet g may by highly inclined to the rest of the system, see errors of in Tab. 2.

These results confirm a hypothesis in the discovery paper. In accord with this work, planetary orbits in the Kepler–11 system should be mutually inclined by no more than a few degrees. It flows from estimating a probability that for a given orientation of the orbits, all six planets transit the star. This reasoning assumes that all inclinations are independent. However, we found that Kepler–11 system is composed of two or three sub-systems, which exhibit small mutual inclinations of orbits. Although a probability that the mutual inclinations between these sub-systems are significant seems a bit larger than for fully independent orbits, it still remains very small. We estimate that a randomly located observer can detect transits of all 6 planets with a probability as small as , for both models I and II.

### 3.3 Models III (eg=0,Ωg=0) and IV (eg=0,Ωi=0)

The results for model III and model IV are given in Tabs. 3 and 4. Most of these results are common for all transit models I to IV. There are some differences regarding a determination of the mass of planet g. In the realm of models III and IV (note that both have fixed and ), only an upper limit of smaller than  Earth masses may be found. The low limits of in model I are likely due to weakly constrained and .

Let us recall that in the bootstrap set derived for model II, we found only two solutions with after  yr. However, this integration time scale is too short to detect weak instability which actually leads to catastrophic disruption of these configurations. It was verified by the direct, long–term integrations. Hence, we did not detect any long–term stable configuration in the bootstrap set of model II. Similarly, stability tests performed for configurations of model III did not reveal any stable models. As compared to model II, fixed seem does not change the general view of the stability of the system.

For model IV we found many stable configurations which confirm stability analysis in Lissauer et al. (2011). They found some stable solutions assuming that the Kepler-11 system is strictly co-planar. We may conclude that a factor of small relative inclinations is more important for maintaining the long term stability than small eccentricities. This will be discussed in more detail further in this work.

We examined a probability that a randomly located observer could detect transits of all planets in the system. This is basically unlikely for model III (), while for model IV a probability of such an event is larger, and we estimate it .

## 4 Discussion on the planet interiors

Figure 10 shows bootstrap diagrams of a few selected pairs of parameters. These results are for model II. The top row is for the semi-major axes and the planetary masses, the radii and mean densities, respectively. The red and green curves mark the data for Uranus and Neptune. The bottom row is for the mass–radius, mass–density and radius–density relations, respectively. Similarly, the red and green filled circles are for Uranus and Neptune. As we concluded above, the orbital solutions in set II are only marginally stable, however, it is a matter of unconstrained orbital angles. Note that a discussion in this Section concerns semi-major axes (known with an excellent precision) as well as planetary masses and radii.

This figure reveals that the most inner four planets in the Kepler-11 system exhibit a clear and curious anti/correlation of masses, radii and densities with the semi-major axes. Masses and radii increase with , while densities decrease. The last panel constructed for shows a weak anti-correlation between the radii and densities, the smaller radius, the larger density.

The determined masses and radii of the planets provide some insight into their chemical composition. We use a simple analytic relation between the radius and the mass of a cold body (Lynden-Bell & O’Dwyer, 2001; Lynden-Bell & Tout, 2001) to estimate the characteristic density of planetary matter. This density can be compared with calculated for a given number of nucleons per number of electrons of a chemical mixture forming the planet, . The value of is a simple mean over the elements in each chemical substance or component. For instance the H-He mixture has for the mass proportions to , and ice or rock has . In this way, we can obtain some insight into likely chemical composition of the planets.

Our results for Kepler-11 are presented in Fig. 11. Black curves with grey areas are for and its uncertainty . Each panel is for one planet of the Kepler-11 system. Data for planets b to planet g are displayed from the left to the right, respectively. The colored curves are for the Solar system, i.e., for Uranus (red), Neptune (green), Jupiter (blue), Saturn (violet) and the Earth (light blue). The density was computed in a wide range of . These values are known relatively well for the Sun companions. Following Helled et al. (2011), for Uranus and Neptune one finds (for the icy model) and (for the rocky model). The density in these particular case is plotted with filled circles. Similarly, for Jupiter and Saturn, may be also estimated (Guillot, 1999). Values of for these particular are marked with circles. Let us note that densities of Jupiter and Saturn are almost identical.

Lynden-Bell & Tout (2001) argue that is the zero–pressure density of the chemical mixture of the planets. Because their model has many simplifications, is usually times larger then . Keeping this in mind, the densities calculated for Kepler-11 planets can be compared with those of the Solar system planets. The value of is best determined for planet e. Its is very close to the Jupiter/Saturn (J/S) value . This suggests, that planet e is built mainly of a H/He mixture with mass proportions of the elements close to with a portion of heavier elements contained in ices or rocks. This makes the planet classified as a smaller “cousin” of Jupiter and Saturn rather than of Neptune and Uranus, as suggested in (Lissauer et al., 2011).

The density of planet b is determined worse than for planet e. It is rather unlikely though that it belongs to the same class as planet e. Parameter is larger than for Jupiter and Saturn, even taking into account a large uncertainty. It is also larger than for Uranus and Neptune (U/N)–like planets but is smaller than for the Earth. We can conclude that planet b is a small planet containing a large percentage of heavy elements in its interior, which is likely larger than in the ice giants. It is reasonable to classify this planet in the super-Earths family, although, it may be also a small Neptune–like planet.

Planet f has the mass similar to planet b. However, its composition is likely between the Jupiter–Saturn and Uranus–Neptune classes. Planet d has likely similar composition as planet f. The best-fit estimate of for planet c is very close to the Uranus/Neptune value. For planet g there is only the upper limit of . However, it is probably close to the Jupiter/Saturn value or, less likely, to the Neptune/Uranus value.

These conclusions are reinforced after inspecting the bottom-left panel of Fig. 10 illustrating the mass–radius diagram for the Kepler-11 system. The mass–radius relation computed for and of the Solar system planets are plotted with different colors. Data are shown for Uranus (red), Neptune (green), Jupiter (blue), Saturn (violet) and Earth (light blue), respectively. Representations of this relation are plotted in the mass-density diagram (the middle panel) and the radius–density diagram (the right panel).

For small masses, the characteristic density is very close to (see Eq. 34 in Lynden-Bell & O’Dwyer, 2001). This derived decay of with the mean distance from the star suggests that the inner planets may contain larger amount of heavy elements than more distant companions. If this correlation can be confirmed, it may provide an observational constraint for the planet formation theory. Allowing for some speculations here, let us note that all Kepler-11 planets exhibit masses in the same range of a few Earth masses. Hence, they likely have formed in a similar way and physical environment (Rogers et al., 2011). Small eccentricities and small relative inclinations suggest that the system evolved orbitally smoothly towards the current state, conserving the ordering of initial distances from the star. The observed relation between and may then indicate the chemical composition and mass density distribution in the primordial protoplanetary disk.

We underline that the results in this section must be considered as preliminary. Due to relatively narrow time–window of the photometric data, masses of the planets are determined with large uncertainties.

## 5 Results of the dynamical analysis

The best–fit solutions gathered with the help of the bootstrap algorithm provide us primary information required to perform extensive study of the dynamical stability of the system. Because the orbits of Kepler-11 super-Earth planets are confined within the mean distance of Mercury in the Solar system, we could expect that the dynamics of this system are extremely complex. Indeed, preliminary integrations demonstrated that the Kepler-11 system is dynamically packed, in accord with a definition and the PPS hypothesis in (Barnes et al., 2008). In spite of apparently ordered configurations with quasi–circular, almost coplanar orbits, and relatively small masses, no long-term stable model I solutions were found. Below, we try to resolve this paradox and we try to detect sources of this seemingly odd and strong instability. To illustrate the structure of the phase space close to the best–fit configurations, we choose a few representative solutions and we construct the MEGNO maps in their vicinity.

### 5.1 Quasi-stable solutions in transit model II

For model II, among initial conditions, we found only 2 configurations exhibiting MEGNO close to after  yr. Parameters of these solutions are listed in Tabs. 5 and 6.

The first stable solution (refereed to as IIa from hereafter, see Tab. 5), has a relatively low mass of planet c Earth masses. Two innermost planets b and c have almost coplanar orbits. The next three planets, d, e and f also form a nearly coplanar sub-system (d-e-f), which is inclined to the first two orbits by large angle . The outermost orbit is inclined even more, by to the inner subsystem of (b-c), and by to the triple–planet subsystem of (d-e-f).

The second stable solution IIb (see Tab. 6) has all masses close to the nominal best–fit values. The mutual inclinations of five inner orbits are close to , while the outermost orbit of planet g is highly inclined to the inner orbits by , similarly to solution IIa. Because the relative inclinations between particular pairs of orbits are large in these best-fit solutions, such systems might be unlikely observed by a randomly located observer. We estimate a probability of such an event as and for solutions IIa and IIb, respectively.

### 5.2 Triple-planet resonances as the main source of instability

Let us now study the vicinity of these particular solutions through the dynamical maps. For each initial condition of the discrete grid with  resolution, we compute over  yr. Figure 12 shows the MEGNO maps for solution IIa. Each panel is for a different pair of planets. The coordinate axes are rescaled mean motions centered at their nominal :

 xi≡ni−ni,0ni,0×104.

The –axes span uncertainties of the semi–major axes , in accord with Tab. 2. The semi-major axes are determined very precisely, hence the interval span a range of to . The rescaled are confined to order of .

The MEGNO is color-coded in the dynamical maps. Blue regions mean regular solutions with , while yellow color is for chaotic (unstable) initial conditions, . The resolution is pixels, total integration time is  yr per pixel, SABA integrator step-size is 0.5 days. A single computation of each map took  hrs of 1200 CPU cores. The integrations of each pixel were performed up to the end time , regardless a value of MEGNO.

Still, although the maps cover tiny regions of the phase space, close to the fixed initial condition, they reveal a sophisticated structure. Because we consider the dynamics in terms of conservative, close to integrable Hamiltonian system, this structure is governed by the resonant motions. A relatively short integration time characteristic periods, equivalent to the orbital period of the outermost planet makes it possible to detect unstable MMRs. They appear as yellow straight bands of different widths and slopes. Inspecting the dynamical maps, we can identify particular resonances.

A condition for the mean motion resonance in the -planet system may be written in the following form:

 N∑i=1pidλid\,t=O[fω,fΩ],orN∑i=1pini=O[fω,fΩ], (3)

where is the mean motion of the -th planet, and are the fundamental frequencies associated with the pericenter arguments and the longitudes of nodes (for all orbits), and are relatively prime integers. The linear relations must obey the d’Alambert rule.

The two-planet MMR takes place when two values of are non-zero. If three coefficients in this linear relation are non-zero, it means that the system exhibits 3–body MMR. In the Kepler-11 system, the fundamental frequencies associated with and are much smaller than (these are the secular frequencies), hence the right-hand sides of Eq. 3 are close to . This makes it possible to skip the secular terms, as the first order approximation, and to identify the MMRs through approximate resonance conditions involving the mean motions only.

To identify the MMRs in the MEGNO maps, we apply a simple method described in (Guzzo, 2005). In the –plane444Let us note that Fig. 12 shows the -planes but may be easily computed from these data., the slope

 αj,k≡dnkdnj

of a particular resonance line determines the ratio of coefficients , i.e.,

 αj,k=−pjpk. (4)

If then the MMR forms a horizontal line, planet is involved in the MMR, while planet is not. If then the MMR forms a vertical line, planet is involved in the resonance, while planet is not. In these cases, other planets may be involved in this particular resonance. If is finite and non–zero then both considered planets are involved in the MMR. To identify this particular resonance, one has to compute slopes corresponding to this resonance in all -planes. It may be not possible, if the map ranges do not cover the resonance band. If is non-zero and finite only for one pair of planets (, ), it means that 2–planet MMR is present. It should be verified whether . Coefficients , of the MMR condition can be computed from the slopes . Similarly, the 3-body MMR takes place, if there exist relatively prime integers , and , such that and are all finite and non–zero. The 3-body MMR may be identified by computing the slope coefficients in appropriate planes of the mean motions. An identification of –body and –body MMRs can be derived as well.

Using this simple MMR identification algorithm, we found most significant MMRs close to solution IIa. The identified 3–body MMRs were labeled at the panels, and listed in Tab. 7. The mean motions of solution IIa permit a few low–order 2–body MMRs in the vicinity of this solution, e.g.,