# Testing a hypothesis of the Octantis planetary system

###### Abstract

We investigate the orbital stability of a putative Jovian planet in a compact binary Octantis reported by Ramm et al. We re-analyzed published radial velocity data in terms of self-consistent Newtonian model and we found stable best-fit solutions that obey observational constraints. They correspond to retrograde orbits, in accord with an earlier hypothesis of Eberle & Cuntz, with apsidal lines anti-aligned with the apses of the binary. The best-fit solutions are confined to tiny stable regions of the phase space. These regions have a structure of the Arnold web formed by overlapping low-order mean motion resonances and their sub-resonances. The presence of a real planet is still questionable, because its formation would be hindered by strong dynamical perturbations. Our numerical study makes use of a new computational Message Passing Interface (MPI) framework Mechanic developed to run massive numerical experiments on CPU clusters.

###### keywords:

-body problem, numerical methods, binaries## 1 Introduction

The Octantis is a single–line spectroscopic binary composed of the Oct A K1 III giant primary ( M) and an unseen red dwarf secondary Oct B, K7-M1 V ( M) separated by au. Colacevich (1935) collected 11 radial velocities (RVs) of the primary, and Alden (1939) determined the astrometric orbit. Recently, the orbital period days and eccentricity were determined by Ramm et al. (2009). Remarkably, they also derived the inclination with an error less than on the basis of Hipparchos astrometry and new 223 precision RVs. Ramm et al. (2009) discovered residual variability of these RVs at a level of and explained this by an S-type (Dvorak, 1984) Jovian planet , Oct Ab , having a minimal mass m in an orbit with a semi-major axis au and an eccentricity .

The putative planet is unusual, because the derived semi-major axis implies its orbit roughly in the middle between massive primary and secondary. A formation of such a planet in the compact binary can be hardly explained (e.g., Kley, 2010). Analytic or general stability criteria formulated for the restricted and general three-body problem like the Hill stability criterion, the results of Holman & Wiegert (1999), or the resonance-overlap criterion by Wisdom (1983) are violated if the planet is assumed to be in a prograde orbit. According with the determination of stability limits in binaries (Holman & Wiegert, 1999), a coplanar, prograde planetary configuration might be stable only up to astrocentric distance au (Ramm et al., 2009). Indeed, Eberle & Cuntz (2010) found that the Keplerian best-fit model in the discovery paper is strongly unstable in just 1–10 binary periods time-scale for the case of prograde motion. This contradicts the planetary hypothesis, although the discovery paper basically excludes non-planetary sources of the observed signal, like stellar spots or pulsations.

Many explanations of the observed RVs variability are still possible, which might resolve the apparent paradox of the strongly unstable system. The residual RVs may be implied by stellar chromospheric activity, a different number of planetary companions, systematic errors of the observations, instrumental instabilities and non-Gaussian uncertainties, such as the red noise (Baluev, 2011), or a different orbital model than the configuration in the discovery paper. Indeed, Morais & Correia (2012) consider a hierarchical triple system with an unseen binary in the place of secondary Oct B. Such an orbital setup forces a precession of the primary orbit around the barycentre of the inner binary and it could mimic the RV variations attributed to the putative planet.

Earlier, Eberle & Cuntz (2010) following the stability studies of the restricted three body problem by Jefferys (1974) found that a retrograde orbit of the planet lies in much wider stable zone than the direct configuration. In a new work, Quarles et al. (2012) investigate this model, extending the set of orbital configurations, through a numerical study in the framework of the elliptic, restricted three-body problem.

In more than a decade, modeling the precision RVs teaches us that the dynamical analysis of the Octantis system done in the discovery paper (Ramm et al., 2009), and following studies (Eberle & Cuntz, 2010; Quarles et al., 2012) do not seem fully consistent with its dynamical character. Because the secondary mass may be almost a half of the primary mass, the relatively wide planetary orbit is strongly perturbed. The perturbation parameter, expressed by the mass ratio of the binary, can be as large as 0.38 (Ramm et al., 2009). In such a case the mutual interactions are important to compute the RV signal (Laughlin & Chambers, 2001). Furthermore, it was assumed that the whole system is co-planar, no matter if the planetary orbit is prograde or retrograde. However, non-planar orbits in compact binaries might appear due to violent post–formation scenarios due to the planet–planet scattering (e.g., Adams & Laughlin, 2003). This may imply highly inclined configurations and their dynamics much complex than in a coplanar case (e.g., Migaszewski & Goździewski, 2011). The dynamical simulations done by Eberle & Cuntz (2010) concern initially aligned orbits. Quarles et al. (2012) study a system with initial planetary eccentricity and two initial longitudes ,, respectively (corresponding to 3 o’clock and 9 o’clock positions in their terminology). While a particular initial orientation of the orbits may be protecting factor for the stability, the initial orbital phases should not be fixed or selected a’priori, if the tested configuration is meant to be consistent with the observations, where we aim to analyze the dynamics of a real system. As a simple example, consider Keplerian RV signal computed from the well known formulae where is the semi-amplitude, is for the pericenter argument of the orbit, is the eccentricity and is the true anomaly. Function changes its sign, if instead of one sets at the initial time. Obviously, a more complex modification of is introduced if other angles are altered, such as orbital longitudes. In general, the initial orbital configuration determines the dynamical evolution and the reflex motion of the primary, which should be consistent with observations.

In this work, we aim to verify and improve the kinematic (Keplerian) model of
the Oct planetary system by searching for the best-fit configurations in
terms of self-consistent dynamical, -body model (Laughlin &
Chambers, 2001, and references
therein). To resolve the fine structure of the
phase space, and to investigate the dynamics and long–term stability of the
planet, we apply the fast-indicator MEGNO (Cincotta et al., 2003) adapted to our
new multi-CPU computing environment Mechanic. In this paper, we
demonstrate a non-trivial application of this software. It was announced by
Słonina et al. (2012).
The Mechanic code is described in
an accompanying work to this paper (Słonina et al., in preparation), and is
freely available on-line (`http://git.ca.umk.pl`

).

This paper is structured as follows. After this introduction, section 2 is devoted to the dynamical analysis of orbital solution proposed in (Eberle & Cuntz, 2010) to possibly global extent. We attempt to visualize complex resonant structures in the vicinity of presumable retrograde planetary configurations. This part has a general character, because we analyze unstable behaviours of S-type planets in compact binaries. In the next section 3, we focus on an re-analysis of the RV observations and on deriving the best-fit parameters of the Oct system in terms of three distinct models of the RV: the Keplerian (kinematic) formulation, Newtonian (dynamic) model, and the Newtonian model with stability constraints (GAMP, Goździewski et al., 2008). In section 4 we demonstrate that stable, retrograde orbits consistent with the RVs data may be found, indeed, but these three RV models lead to significantly different orbital architectures. Conclusions are given in section 5.

## 2 The dynamics of the Octantis system

To resolve the paradox of the Octantis system instability, Eberle &
Cuntz (2010)
searched for stable orbits at a grid of primary mass ( and )
M, three mass ratios ( = , ), and 30–step
initial distance ratio between the planet and the secondary, au. The equations
of motion were integrated for years ( 350 binary periods
) for prograde orbits, and for years ( 3500
) for retrograde orbits. The initial orbits were aligned
with both secondaries fixed in their apoastrons with respect to the primary.
These integrations confirmed the theoretical stability limit au for prograde orbits in (Holman &
Wiegert, 1999). For the retrograde
case, the stability limit was found much larger, indeed,
consistent with the formal best-fit error in
(Ramm et al., 2009)^{1}^{1}1
We found some ambiguity in the calculation of
by Eberle &
Cuntz (2010), who state that in the
elliptical case, denotes the ratio of the initial distance of the
planet from the primary () relative to the semi-major axis of the
binary components (). Following this, for their retrograde
configuration stable for 10 Myr, we may compute
, where ,
au, and , hence au. It might be also au, if the initial apoastron distance of the planet was
approximated as , or, if then au,
which is most likely approximation used in the study of Eberle &
Cuntz (2010).
.

### 2.1 Stability of a system with initially aligned apsides

To extend the study of Eberle &
Cuntz (2010), we applied the
dynamical maps technique. It helps to illustrate the global structure of the
phase space, as well as to identify possible resonances as the origin of
strong instability observed in this system.
We use MEGNO which is a measure of the maximal Lyapunov exponent
(Cincotta et al., 2003) as the stability indicator.
Our serial-CPU software to integrate the equations of
motions and the variational equations of the planetary problem was
encapsulated in the Mechanic MPI module called the CSM^{2}^{2}2
This code is available upon request.. A strongly
interacting planetary system with collisional orbits requires an appropriate
integrator. We applied the Bulirsh-Stoer-Gragg scheme implemented as
reliable and well tested ODEX routine (Hairer
et al., 1993). It provides very
good accuracy and performance (see also, Chambers, 1999). The
total integration time of a single initial condition was not less than
periods of the binary, and in some cases more than
periods. This time is typically one–two orders of magnitude longer than the
time-span of direct numerical integrations in Eberle &
Cuntz (2010), and provides at
least - longer estimate of the Lagrange stability, in accord
corrwith features of MEGNO (Cincotta et al., 2003). The time scale of MEGNO integrations is
long enough to detect most significant mean motion resonances (Goździewski et al., 2008).

In the first numerical experiment, we set the primary mass to M
which is equivalent to . In accord with Eberle &
Cuntz (2010), the planet
and the secondary are placed in the apoastrons of their orbits, hence the
initial mean anomalies , of
the planet and the secondary, respectively. Arguments of periastrons and the
nodal longitudes were set to . In the following numerical
experiments, we varied the inclination of the planet (equivalent to the
relative inclination of the orbits) between (direct
coplanar configuration) and (retrograde coplanar
configuration). This setup, which we consider as aligned orbits,
reproduces the initial configurations in Eberle &
Cuntz (2010), i.e., the 9 o’clock
position in their terminology^{3}^{3}3In the literature, the terms of aligned and antialigned periapses (or orbits, apsidal lines) are well
defined only for direct orbits. Here, for coplanar, direct and retrograde
orbits these terms are defined through the difference of the longitudes of
periapses, equal to and , respectively..
We computed
MEGNO maps in high-resolution with a grid of pixels
that were integrated for
25,000 periods of the binary. The results are illustrated in
Fig. 1, in the semi–major axis — relative
inclination –plane for the two secondary masses.
The right map is for an
artificial brown dwarf secondary, to illustrate, through a comparison, the
strength of the perturbation in the nominal Oct system (the left panel in
Fig. 1). For prograde configurations, the stability limit
au is slightly smaller than its value found in the
previous papers. For retrograde orbits, it is almost doubled and reaches
au. A putative planet could be located at a border of
this zone, filled with unstable low-order MMRs. A stable retrograde orbit
found by Eberle &
Cuntz (2010) lies in a region spanned by strong 4:1, 3:1, 5:2 and
2:1 MMRs. A close-up of this area covering the 5:2 MMR,
slightly beyond the semi-major axis of the planet determined
in (Ramm et al., 2009), is shown in
Fig. 2. Overlapping of the MMRs and their sub-resonances creates
a global chaotic zone. This effect is identified as the main source of
instability of planets in compact binaries (Mudryk &
Wu, 2006).

The stability limit depends strongly on the relative inclination and the mass of the secondary. Configurations with intermediate relative inclinations, roughly between and , and particularly solutions close to polar orbits () are very chaotic. They are associated with the Kozai resonance (Kozai, 1962), see also section 2.3. For non-coplanar orbits, the border of stable zone may be much closer to the primary than 0.4 au. Even for a hypothetical, low–mass secondary 14 m (Fig. 1, right panel), polar orbits are very unstable. In this case, stable prograde solutions may be found up to au. We also note that individual unstable MMRs are much wider for prograde orbits than for retrograde configurations. A theoretical explanation of this phenomenon is given in a recent paper by Morais & Giuppone (2012). The stability limits are governed by the phase-space topology of retrograde and prograde MMRs. At the / mean motion ratio, the prograde resonance is of order while the retrograde resonance is of order , hence the resonance order is higher. Although this result is derived for the circular, restricted three body problem, it should be valid also for small–eccentricity binaries. This is confirmed by Fig. 1 showing different widths of prograde and retrograde MMRs for .

### 2.2 Developing the instability with the secondary mass

An identification of the origin of instabilities for large mass ratio of the binary is not straightforward. Due to overlapping MMRs and their sub-resonances, a complex pattern of stable and unstable motions appears (compare panels in Fig. 1). To study how the unstable zones are evolved with increasing perturbation, we performed a numerical experiment in which the secondary mass was gradually increased from an artificial value of 14 m (a brown dwarf, see the right panel in Fig. 1), up to the nominal value of 532 m (the left panel of Fig. 1). For the small mass secondary, the MMRs can be easily identified, and we can follow, how they expand when the perturbation of the secondary grows. In this way, we may also refer to a low-dimensional dynamical system given through model Hamiltonian in Froeschlé et al. (2000). The dynamics of the S-type planet in the binary should reveal qualitatively similar features because this is also governed by a perturbed Hamiltonian of the form , where is the integrable (Keplerian) part, and is the perturbation due to mutual interactions of the planet and the binary. The perturbation parameter can be expressed through where is the mass of the secondary companion Oct B and is the mass of the primary Oct A.

The results are illustrated through a sequence of dynamical maps in Fig. 3. These MEGNO maps were computed for periods of the binary. At each panel, we labeled the lowest-order MMRs and the mass of the secondary (perturber). Each Mechanic run occupied 256 CPUs, and computations took, depending on initial conditions (an extent of the chaotic zone) up to a few hours of CPU time. Already for the smallest mass ratio, corresponding to brown dwarf mass perturber, a large part of the phase space is strongly chaotic (the top-left panel of Fig. 3, also the right panel in Fig. 1). In the relevant zone of retrograde motions, a pattern of narrow MMRs appears. A retrograde orbit might have much more stable space to explore than a direct orbit. In this regime, the probability of a stable orbit would be close to 1. When the perturber mass is increased, the widths of MMRs quickly expand, and for m (roughly, the left mass limit of the secondary) one may observe that strong MMRs are overlapping which emerges a zone of global chaos, roughly beyond au. For the nominal mass of the perturber (see the bottom right panel in Fig. 3), there are only narrow areas of stable motions.

A close-up of that region in Fig. 2 demonstrates a view of the phase space of the model Hamiltonian for large perturbation parameter (Froeschlé et al., 2000), see also (Słonina et al., 2012). The fine structures seen in this map represent the Arnold web in the three body problem. They appear due to overlapping MMRs and their sub-resonances, strongly split by the perturbation. This feature is expected because it was found in the Outer Solar system (Guzzo, 2005) and, very recently, in the Kepler-11 extrasolar system of six planets (Migaszewski et al., 2012). Likely, it has not been shown in the three body problem with so much details. The right panel of Fig. 2 shows the MEGNO map in the same area, but with lower resolution of pixels, computed for periods of the binary. We see that a fine resolution pixels, as well as longer integration time P are crucial to discover the fine details of the phase space. ( See section 4 for more examples of such structures in the neighborhood of the Newtonian best-fit model of the Octantis system, consistent with available observations).

The dynamical stability of planetary orbits in systems exhibiting strong mutual interactions can be influenced even by very small changes of the initial conditions due to the presence of unstable resonances and their overlapping. If the mutual perturbations are small enough, and chaotic motions appear on a regular net, they may be practically stable over very long times (Guzzo, 2005). Such a state of the system is called the Nekhoroshev regime. Otherwise, if the chaotic motions do not constitute a regular web, and most of orbits form a global chaotic zone, the stability of the system is influenced by a strong chaotic diffusion. This regime is related to the resonance overlapping, and is called the Chirikov regime (Froeschlé et al., 2000; Guzzo, 2005).

The Arnold web emerging due to three–body and four–body MMRs was investigated by Guzzo (2005, 2006) in the Outer Solar system. He mapped the phase space in the semi-major axes planes of Jupiter, Saturn, Uranus and Neptune. Here, we illustrate this feature in the three–body problem in the semi-major axis — inclination plane, corresponding to a different choice of the canonical actions, and appearing due to two–body MMRs. The Arnold web creates an intermediate zone between the ordered and strongly chaotic motions, which spans relatively wide range of the inclination, extending for . In this intermediate region, the phase-space motions which are chaotic, may persist over long periods of time. The results of this experiment also show that the stability limits for co-planar binaries derived numerically by Holman & Wiegert (1999) may be strongly affected by the inclination. In our case of moderate eccentricity of the planetary orbit, these limits are roughly valid for and . A determination of stability limits in systems with non-zero relative inclinations needs additional extensive numerical study recalling that our computations concern a particular initial configuration of the orbits. We postpone this subject to a new work.

### 2.3 The secular dynamics of the planet

The origin of a wide strip of unstable motions around polar orbits can be explained with the help of the secular octupole–level theory (e.g., Migaszewski & Goździewski, 2011, and references therein). For a direct comparison with the MEGNO integrations, we focus again on the initial orbital setup in (Eberle & Cuntz, 2010). For small au and less, the expansion parameter where the system may be considered as non-resonant. This permit us to study the long–term dynamics through the secular theory. The range of au is investigated to compare both relevant intervals of direct and retrograde motion although becomes as large as at the right border of this interval and the assumptions of the secular theory may be questionable.

The results are illustrated in Fig. 4 in the form of dynamical maps in the )–plane. The left column is for the nominal mass of Oct B, and the right column is for an artificial, low-mass secondary of m. The top row illustrates the maximal eccentricity attained by the planet over the secular orbital period. Both perturbers are massive enough to force extreme eccentricity . This must imply collisions with the primaries. Indeed, in the middle panels showing the maximum apocenter distance of the planet, the filled dots mark initial conditions that led to a collision with the star. These two zones closely coincide. If the initial semi-major axis is increased, the minimum distance of the planet to the secondary becomes almost always close to 0 au (see the bottom row in Fig. 4). Filled dots in these panels mark regions, in which angle librates.

This experiment reveals that a continuous unstable zone between predicted by the secular model closely coincides with unstable region around the polar orbits exhibited in the full three-body system. It means that this unstable region appear due to the Kozai resonance which is able to force large eccentricities in the secular time-scale. More details and a concise description of the Kozai resonance in the context of planetary dynamics may be found in a recent review (Beaugé et al., 2012, p. 1065–1068, and references therein). Hence the primary source of instability can be identified with a secular effect rather than with short-term chaotic dynamics associated with the MMRs. Outside unstable zones of the Kozai resonance, the primary source of instability are the MMRs which appear as discrete areas of chaos for small value of the perturbation parameters, and a zone of global chaos when the perturbation becomes sufficiently strong and the short-term MMRs may govern the global dynamics. Let us note, that although the global picture of the secular system is illustrated here for the aligned configuration, it remains qualitatively the same for the anti-aligned orbits.

## 3 Re-analysis of the RV data

It is well known that the RV measurements errors originate from the formal errors of the Doppler method, systematic uncertainties due to instrumental instabilities, as well as from internal chromospheric variability of the parent star. The later error is difficult to quantify particularly for young or active stars. In the recent literature, it is common to derive astrophysical estimates of the jitter (Wright, 2005). Joint error of the measurements are obtained by adding both uncertainties in quadrature for each data point. Recently, Baluev (2009) has shown how to estimate as an additional free parameter of the fit model. Here, we follow the earlier, somehow more traditional approach. In this work, we computed the best-fit models after rescaling the measurement errors in (Ramm et al., 2009) by two values of . For the first class of models, we choose , which is relatively small value, exhibited by quiet, non-active stars, and likely too conservative in our case. The second family of solutions was obtained for . This value of jitter uncertainty seems better justified recalling large residuals of the Keplerian model in the discovery paper of a similar magnitude ( ). Using these RV data, we optimized different 2-planet models. The astrometric data gathered by Hipparchos were used indirectly, by fixing the nodal line and inclination of the binary orbit close to the combined astrometric and radial velocity solution in (Ramm et al., 2009).

### 3.1 Keplerian (kinematic) model of the RV

In the first fitting experiment, we search for the best-fit Keplerian, non-interacting orbits. Hence, we basically follow the discovery paper. The best-fit configurations of Oct system were found after extensive search performed with the hybrid optimization algorithm (see Goździewski et al., 2008, and references therein). It consists on two stages: a quasi-global Genetic Algorithm which provides many reasonably precise models, and the Levenberg-Marquardt (LM) method which is used to refine these solutions.

The results are illustrated in Fig. 5, in the form of parameter scans illustrated in the orbital period – eccentricity plane of the planet. Curiously, for both selected jitter uncertainties, we found a well defined secondary minimum for days and , besides the best-fit model reported by Ramm et al. (2009) with days and . A choice of the jitter uncertainty may lead to significant changes of the shape and alters formal errors of the best-fit parameters. This effect is important for fixing parameter bounds in dynamical studies made after the best-fit models are found. Moreover, because the secondary best-fit configuration has the planetary period comparable with the binary period, its interpretation in terms of the osculating elements becomes even more questionable than in the case of the primary minimum of , as underlined in the introduction. The secondary solution might correspond to 1:1 MMR which obviously would require a more general -body model of the RV which takes into account all mutual interactions between the bodies in the system (Laughlin & Chambers, 2001).

### 3.2 The self-consistent, Newtonian model of the RV

As we have shown above, the Keplerian model of Oct data is non-unique, permitting two different solutions. It is also well known that it is independent on orbital inclinations and true masses. Due to the massive secondary, the RV-model including mutual interactions between all system components might constrain these parameters, in spite of a narrow observational window, like in other interacting extrasolar planetary systems (Laughlin & Chambers, 2001). Hence, in the second fitting experiment, to compute the RVs of the primary, we integrated the -body equations of motion. In this model, the primary mass , its semi–major axis , eccentricity of the binary, pericenter argument , inclination and the initial mean anomaly of the binary were allowed to vary within the formal error bounds, in accord with the combined RV and astrometric solution in Ramm et al. (2009). The nodal longitude of the binary was fixed, . The formal error of this element is irrelevant for the dynamics of the system, because selecting the nodal line of the binary corresponds to a choice of the –axis of the reference frame. The planet mass , its inclination , the nodal longitude and the remaining orbital elements with the RV offset were added to the set of 14 free parameters of the fit model. The uncertainties of the measurements in (Ramm et al., 2009) were rescaled in quadrature, as in the Keplerian model, to account for the intrinsic stellar RV variability.

Similar to the Keplerian model, we explored the parameter space with the hybrid algorithm. After an extensive search, we selected the best-fit -body configuration. Then we reconstructed the shape of around this nominal solution in 2-dimensional planes of various orbital elements. Its parameters derived with are given in Table 1 (Fit I). The RV data in Ramm et al. (2009) consist of many “clumps” of measurements done during one night, which cover less than 0.1 day. Such data points were binned, and the resulting data set consists of 91 measurements. This step was helpful to improve the CPU-efficiency, due to necessity of computing numerical derivatives in the LM algorithm. We examined that -scans can be reproduced with the full, original data set. Obtaining a smooth shapes of for non-binned data is difficult, due to numerical errors introduced by the numerical derivatives.

The –plane (Fig. 6, the top-lef panel) makes it possible to compare the topology of in both RV-models. Clearly, the best-fit Newtonian configuration is qualitatively different from the Keplerian models found in the previous section which confirms our earlier predictions. The planet eccentricity may be as large as . The low estimate of implies two, statistically distinct close minima of , for au, and better model for au. These minima are resolved at the formal –error level of . If the jitter estimate is larger, likely much more realistic, only one solution persist (see Fig. 7, the top left panel). In this case, the formal error of spans almost au.

Even more interesting results are illustrated in the next panels of Figs. 6 and 7 showing scans of in the semi-major axis–inclination, —longitude of the node and —the pericenter argument planes, respectively. scans computed for reveal two minima seen in Fig. 6. Such scans for (Fig. 7) illustrate that the global minimum around au, is well defined in all planes of angular elements. This means that even a narrow observational window, covering roughly 2 planetary periods, makes it possible to estimate all angles of the orbit due to strong mutual interactions in the system. This experiment is a convincing illustration that the Keplerian, kinematic model of Oct is questionable because it looses information of the node and inclination of the orbit. A proper estimate of the jitter uncertainty is vital to derive true errors of the best-fit model. We found that the best-fit model (Fit I in Table 1) tends towards the lower bound of inclination (69), with simultaneous increase of the secondary mass. This effect appears because the inclination of the binary cannot be constrained by the RV data alone. Likely, dynamically consistent –body model of the Oct system should also incorporate the Hipparchos astrometry, which would be helpful to fix the orbit of the binary fully consistent with all data and the three-body model. We postpone this analysis to another work.

### 3.3 Newtonian model with stability constraints

Unfortunately, the best-fit Newtonian model (Fit I in Table 1) found through the hybrid search leads to strongly unstable system. Due to relatively large parameter error, stable configurations might be still possible in its neighborhood. Hence, we conducted an additional, extensive search for stable models imposing stability constraints (GAMP algorithm, Goździewski et al., 2008). In this approach, unstable models are “penalized” by an artificial, large value of (or an rms) and are excluded during the optimization. The results are shown in Fig. 6. Stable fits exhibit systematically larger than the best–fit model in Table 1. To illustrate the statistics of such configurations, we over-plotted solutions with an rms (white, filled circles) at the top of selected parameter scans of (Fig. 6). Let us note that stable models have the MEGNO signature for . This relatively short integration period requires an acceptable CPU-overhead from the hybrid algorithm, still providing at least 10-100 longer Lagrange stability time. Such solutions should survive for – or even longer integration time (Cincotta et al., 2003). The statistics in Fig. 6 reveals that the formal best-fit configurations in terms of “pure” -body model are systematically shifted by more than from the best-fit GAMP (stable) models. This regards all orbital parameters, and in particular, the eccentricity and the argument of pericenter. This large discrepancy may indicate that the 2-planet model is in fact not proper, and should be rejected at all or modified. Curiously, we found stable solutions only as retrograde orbits, actually in accord with the original hypothesis in (Eberle & Cuntz, 2010) and later investigated by Quarles et al. (2012). This is illustrated in Fig. 8 showing that differences of the longitudes of periastrons and of the nodal longitudes are well bounded around , resulting in apsidal lines of the orbits which are significantly anti-aligned (the left panel). Large relative inclinations of these stable solutions (the right panel) mean retrograde orbits in the range of around 1.2 au.

To compare the quality of these two families of models in terms of an rms, the RV residuals of two selected best–fit models are displayed in Fig. 9: the first one derived without stability constraints (Fit I in Tab. 1, the left panel), and an low rms stable model ( the right panel, Fit II in Table 1). In both cases, the residuals exhibit a large magnitude of – , comparable with the semi-amplitude of the RV signal due to the putative planet itself. In this sense, stable solutions might be still acceptable, although their rms is significantly larger from the rms of the formal best-fit model configuration. A large variability of the RV residuals cannot be well explained by the 1-planet model. This may indicate strong stellar activity mimicking planetary signal.

Fit I (-body, unstable) | Fit II (GAMP, stable) | Fit III (GAMP, stable) | ||||

Parameter | planet OctAb | Oct B | planet OctAb | Oct B | planet OctAb | Oct B |

mass [m] | 3.54 0.07 | 572 15 | 1.92596 | 560.62366 | 2.135457 | 559.2707 |

semi-major axis [AU] | 1.41 0.02 | 2.533 0.005 | 1.16365 | 2.52813 | 1.160281 | 2.526478 |

eccentricity | 0.40 0.01 | 0.2331 0.0006 | 0.13139 | 0.23881 | 0.247587 | 0.238713 |

inclination [deg] | 84.9 2.2 | 69.0 3.3 | 110.21102 | 71.28090 | 102.4152 | 71.681 |

longitude of the node [deg] | 289.2 4.2 | 87.0 (fixed) | 262.29811 | 87.0 (fixed) | 260.0224 | 87.0 (fixed) |

argument of pericenter [deg] | 340.6 2.6 | 75.0 0.1 | 127.26299 | 74.59137 | 88.86501 | 74.85992 |

mean anomaly [deg] | 289.2 4.2 | 339.5 0.l | 133.16327 | 339.762973 | 165.3685 | 339.5231 |

2.79 | 3.47 | 3.44 | ||||

RV offset [m s] | -6417 1 | -6473.682 | -6475.317 | |||

rms [m s] | 19.9 | 25.2 | 24.8 |

### 3.4 Stability of the best-fit GAMP configurations

To illustrate the dynamical neighborhood of the best–fit GAMP models, we selected two examples of low rms configurations. Their osculating elements at the epoch of the first observation in (Ramm et al., 2009) are given in Table 1 (Fits II and III). Many digits are quoted, to reproduce elements of these fit possibly exactly, to make it possible to reproduce their stable evolution due to complex and chaotic neighborhoods. The formal errors may be estimated graphically in Fig. 6 or interpreted as uncertainties of the formal Fit I. The primary mass is fixed as equal to 1.4 M. An rms of this solution is .

The top-left panel of Fig. 10 illustrates the MEGNO dynamical map in the –plane computed for Fit II (see Tab.1). This map shows a global view of the phase space. Note that inclination is absolute w.r.t. the astrometric frame, in which the binary inclination is bounded around . This figure reveals that the phase space is mostly strongly chaotic. The best–fit configuration is found in a narrow stability island close to the 7:2 MMR. For a reference, other significant MMRs are labeled. A close-up of this map shown in the top–right panel of Fig. 10 reveals an extremely complex structure of this regular island. It spans a tiny region covering au and . In this small island, the phase space has again a sophisticated, fractal-like structure of the Arnold web.

The integration time used to derive the global dynamical map in the left panel of
Fig. 10 is . Subsequent close-ups shown in
Fig. 10 were integrated for up to
per each point, depending on the magnification factor. To
detect weak MMRs and their sub-resonances, such long integration times are
indispensable. The resolution of these maps is pixels.
Calculations were performed with the help of the Mechanic installed on the
UV SGI supercomputer chimera of the Poznań Supercomputing Centre
(Poland)^{4}^{4}4Supercomputer chimera has 2048 CPU cores, and runs
spanning up to 24 hrs occupied all available hardware resources. We did not
observed any bottlenecks in the master–worker inter-communication. This nice
performance of the Mechanic code was expected thanks to a small overhead of
the MPI data broadcasts. These computations validated the software and its
ability of stable long-runs on a few thousands of CPU cores..

We computed a similar set of dynamical maps for another stable model, with larger eccentricity (Fit III in Table 1). The results are illustrated in Fig. 11. The top left panel illustrates the dynamical map in the ()–plane, showing a narrow stable island of this solution, the rest of plots are for the ()–plane and its close-ups.

To reveal the Arnold web in more detail, we computed a number of close-ups shown in panels of Figs. 10–12. These maps reveal many weak sub-resonances, which form a regular net present at any magnification level. These structures are particularly well seen in the bottom-right panel of Fig. 12. They closely resemble the Arnold web of the model Hamiltonian (Froeschlé et al., 2000). Likely, it has been never illustrated with such a sharpness and details. We also tested other solutions from the GAMP–derived set. The results are similar. The best–fit models which obey reasonably well the observational constrains are confined to only narrow stability islands.

## 4 Conclusions

In this work, we aim to improve the original Keplerian models of the Oct system. We re-analysed the RV observations in (Ramm et al., 2009). The results confirm our predictions on a significance of the mutual interactions between the secondary and the putative planet. These interactions are so strong that – if the planet is real – it might be possible to derive the orientation of the orbit with respect to the binary, and its true mass, in spite of a narrow observational window spanning less than 2 orbital periods.

We confirm that stable retrograde solutions in the Oct planetary system are possible as proposed in (Eberle & Cuntz, 2010; Quarles et al., 2012). However, the observational constrains imply such orbits with apsidal lines initially anti-aligned with the apses of the binary. The long-term stable solutions may be found only in relatively small islands of regular, quasi-periodic motions exhibiting a complex structure of the Arnold web. Stable models must be significantly shifted in the parameter space from the -body and Keplerian best-fit models. It means that the 2-planet, S-type configuration is most likely inadequate to fit the available data set. It remains highly uncertain how a massive, Jovian planet could be trapped in such tiny stable regions, or how it might form in globally unstable dynamical environment. The observational RV signal is very noisy, and the best-fit models reveal a large scatter of residuals having amplitude comparable with the RV signal itself. These arguments make the hypothesis of Jovian planet in orbit around Oct A questionable. Actually, this in accord with the doubts raised in the discovery paper (Ramm et al., 2009). New observations are required to confirm or withdraw that explanation of the observed RV residual signal, and to search for alternate dynamical models of the detected variability.

An efficient analysis of the observations of extrasolar planetary systems involves not only different observational techniques and data sources, but also orbital optimization algorithms, combined with analytical and numerical stability studies. In this work, we propose and apply our new MPI-based task management system Mechanic to perform high-resolution dynamical mapping of the phase space. This framework is devoted to numerical analysis of large sets of initial conditions, like the long-term integrations of the equations of motion. We work on other applications, like quasi-global optimization with the Genetic Algorithms of different observational models of the extrasolar planetary systems. The Mechanic code with detailed technical documentation and sample modules is freely available at the project website (Słonina & Goździewski, 2012).

## 5 Acknowledgments

We thank the anonymous referees for their reviews that improved this paper. This work is supported by the Polish Ministry of Science and Higher Education through grant N/N203/402739. Computations were conducted thanks to the POWIEW project of the European Regional Development Fund in Innovative Economy Programme POIG.02.03.00-00-018/08.

## References

- Adams & Laughlin (2003) Adams F. C., Laughlin G., 2003, Icarus, 163, 290
- Alden (1939) Alden H. L., 1939, AJ, 48, 81
- Baluev (2009) Baluev R. V., 2009, MNRAS, 393, 969
- Baluev (2011) Baluev R. V., 2011, Celestial Mechanics and Dynamical Astronomy, 111, 235
- Beaugé et al. (2012) Beaugé C., Ferraz-Mello S., Michtchenko T. A., 2012, Research in Astronomy and Astrophysics, 12, 1044
- Chambers (1999) Chambers J. E., 1999, MNRAS, 304, 793
- Cincotta et al. (2003) Cincotta P. M., Giordano C. M., Simó C., 2003, Physica D Nonlinear Phenomena, 182, 151
- Colacevich (1935) Colacevich A., 1935, PASP, 47, 87
- Dvorak (1984) Dvorak R., 1984, Celestial Mechanics, 34, 369
- Eberle & Cuntz (2010) Eberle J., Cuntz M., 2010, ApJL, 721, L168
- Froeschlé et al. (2000) Froeschlé C., Guzzo M., Lega E., 2000, Science, 289, 2108
- Goździewski et al. (2008) Goździewski K., Breiter S., Borczyk W., 2008, MNRAS, 383, 989
- Guzzo (2005) Guzzo M., 2005, Icarus, 174, 273
- Guzzo (2006) Guzzo M., 2006, Icarus, 181, 475
- Hairer et al. (1993) Hairer E., Norsett S. P., Wanner G., 1993, Solving Ordinary Differential Equations I. Nonstiff Problems. Springer Verlag
- Holman & Wiegert (1999) Holman M. J., Wiegert P. A., 1999, AJ, 117, 621
- Jefferys (1974) Jefferys W. H., 1974, AJ, 79, 710
- Kley (2010) Kley W., 2010, in K. Go zdziewski, A. Niedzielski, & J. Schneider ed., EAS Publications Series Vol. 42 of EAS Publications Series, The formation of massive planets in binary star systems. pp 227–238
- Kozai (1962) Kozai Y., 1962, AJ, 67, 591
- Laughlin & Chambers (2001) Laughlin G., Chambers J. E., 2001, ApJL, 551, L109
- Migaszewski & Goździewski (2011) Migaszewski C., Goździewski K., 2011, MNRAS, 411, 565
- Migaszewski et al. (2012) Migaszewski C., Słonina M., Goździewski K., 2012, MNRAS, 427, 770
- Morais & Correia (2012) Morais M. H. M., Correia A. C. M., 2012, MNRAS, 419, 3447
- Morais & Giuppone (2012) Morais M. H. M., Giuppone C. A., 2012, MNRAS, 424, 52
- Mudryk & Wu (2006) Mudryk L. R., Wu Y., 2006, ApJ, 639, 423
- Quarles et al. (2012) Quarles B., Cuntz M., Musielak Z. E., 2012, MNRAS, 421, 2930
- Ramm et al. (2009) Ramm D. J., Pourbaix D., Hearnshaw J. B., Komonjinda S., 2009, MNRAS, 394, 1695
- Słonina & Goździewski (2012) Słonina M., Goździewski K., , 2012, The Mechanic User Guide, http://git.astri.umk.pl
- Słonina et al. (2012) Słonina M., Goździewski K., Migaszewski C., 2012, in Arenou F., Hestroffer D., eds, Proceedings of the workshop ”Orbital Couples: Pas de Deux in the Solar System and the Milky Way”. Held at the Observatoire de Paris, 10-12 October 2011. Editors: F. Arenou, D. Hestroffer. ISBN 2-910015-64-5, p. 125-129 Mechanic: a new numerical MPI framework for the dynamical astronomy. pp 125–129
- Szebehely (1967) Szebehely V., 1967, Theory of orbits. The restricted problem of three bodies
- Wisdom (1983) Wisdom J., 1983, Icarus, 56, 51
- Wright (2005) Wright J. T., 2005, PASP, 117, 657