## Abstract

The three-planet extrasolar system of HD 181433 has been detected with HARPS. The best-fitting solution, announced by the discovery team, describes a highly unstable, self-disrupting configuration. In fact, a narrow observational window, only partially covering the longest orbital period, can lead to solutions representing unrealistic scenarios. Taking into account the *dynamical stability* as an additional observable while interpreting the radial velocity (RV) data, we can analyse the phase space in a neighbourhood of the statistically best fit and derive dynamically stable configurations that reproduce the observed RV signal. Our Newtonian stable best-fitting model is capable of surviving for at least 250 Myr. The two giant companions are found to be locked in the 5:2 mean motion resonance (MMR) as Jupiter and Saturn in the Solar system. This mechanism does not allow close encounters even in case of highly eccentric orbits. Moreover, planets *c* and *d* are located in regions spanned by many other strong low-order MMRs. We study the dynamics of some plausible scenarios, and we illustrate the behaviours caused by secular apsidal resonances and MMRs. Furthermore, we find a terrestrial planet in the habitable zone of HD 181433 can retain stability. Apart from filling an empty gap in the system, this body could offer a harbour for life indeed. Additional measurements are necessary in order to investigate this hypothesis and can confirm the predictions outlined in the paper.

## 1 INTRODUCTION

Nowadays, more than 50 multiplanet systems are known to harbour two to six planets.1 The recent years have seen a proliferation of multiple-planet systems thanks to the increase in both instrument precisions and duration of several planet search programmes. This has allowed the detection of longer periods planetary signatures, as well as planetary signatures with lower amplitude.

An increasing sample of found planets improves our knowledge of their distribution in the mass–period diagram and allows comparison with theoretical predictions (e.g. Wright et al. 2009). In fact, dynamical analysis of planetary systems can both precisely determine their orbital architectures and constrain their evolutionary histories providing a test bed for planetary formation and evolution theories. To make these investigations, we need measured orbital parameters as accurate as possible. Unfortunately, the accuracy of such measurements are limited due to uncertainties and degeneracies inherent to the radial velocity (RV) discovery technique. Nevertheless, the RV method is the most efficient technique for detecting extrasolar planets with more than 90 per cent of all currently known planets being either detected or characterized using this method.

As the time baseline becomes larger, it is possible to distinguish a trend in the RV signals due to long-period outer companions that have not completed a single orbit yet. In this case, it can happen either that the profile of χ^{2} is very smooth or that it does not have a well definite minimum, as a result the confidence levels may comprise large intervals of the parameters fitted.

The RV signal does not provide any direct information on the real masses and the orbital inclinations of the planets; even if we make the assumption the system is coplanar and seen edge-on, an *N*-planet configuration is described by some 5*N*+1 parameters. The a priori unspecified number of planets, narrow observational windows, stellar jitter and weakly constrained orbital parameters can lead to not unique solutions and (quite often) to best fits representing unrealistic scenarios. In fact, these solutions can present well-constrained minimum masses but also a poorly constrained eccentricity for the outermost planet; i.e. in the statistically optimal best-fitting solutions, the eccentricities can be large and can rapidly (on the time-scale of thousands of years) bring to catastrophic collisional instabilities (Goździewski, Migaszewski & Musielinski 2008).

According to the Copernican assumption that we are not observing the Universe at a privileged time, the detection of a rapidly unstable system during a few years of RV observations is not likely, then a fit which corresponds to a quickly unstable configuration is also doubtful. Therefore, we can aim to put limits on the masses and orbital elements of the planets by investigating and finding the plausible and stable solutions. In this logic, the *dynamical stability* is an additional observable that must be taken into consideration when interpreting the RV data. It turns into a discriminating element especially when the longest orbital period is only partially covered (Murray & Holman 2001; Goździewski, Konacki & Maciejewski 2003).

Often, a stability criterion is applied once the best-fitting unstable solution is found, subsequently the orbital elements are tuned to get a stable configuration. However, this approach does not necessarily provide stable fits that are simultaneously optimal in terms of χ^{2} or rms. Most of the times, with the term *stable* we indicate a configuration which does not disrupt or change qualitatively during a period of time of the order of million years. The literature is plenty of studies which take into account stability when modelling the RV data. Here we just mention the works of (1) Veras & Ford (2010) reporting on the system stability, secular evolution and the extent of the resonant interactions for five dynamically active multiplanet systems; (2) Wright et al. (2009) which derive updated orbital parameters for a number of systems considering mutual interactions between planets; (3) Goździewski et al. (2003), Goździewski, Konacki & Maciejewski (2006) and Goździewski et al. (2008) which directly eliminate unstable solutions during the fitting procedure using gamp and (4) Correia et al. (2010) which give constrains on the inclination with respect to the line of sight for some of the planets in the GJ 876 system. We also should remind the study about the directly imaged system HR 8799, where the difficulties in finding regions of stable motion may indicate the system is undergoing a phase of planet–planet scattering (Goździewski & Migaszewski 2009).

Thanks to powerful instrumentations like the HARPS spectrograph at La Silla Observatory, the RV accuracy is increasing rapidly breaking the barrier of 1 m s^{−1} (Pepe et al. 2003). This has permitted the detection of weaker signals consenting the discovery of some of the lowest-mass planets identified such as GJ 876 d (Rivera et al. 2005), HD 40307 b (Mayor et al. 2009a), 61 Vir b (Vogt et al. 2010) and GJ 581 e (Mayor et al. 2009b).

The planetary system of HD 181433 has been discovered with HARPS. It has been reported to contain three planets: two Jupiter-class planets and a super-Earth of 7.5 m_{⊕} (Bouchy et al. 2009).

Inspired by the peculiar properties of the system, which includes two giant planets and one rocky planet all in high eccentric orbits, we aimed to study the past and future evolution of the system. Unfortunately, the published configuration is unstable.2 The model in which the initial eccentricities of the planets are reduced by 1σ quickly leads to disruption too. The fate does not change when we assume, in addition, a mutual orbital inclination of 20° between planets *c* and *d*.

These attempts evidence the necessity of doing an analysis from scratch in order to get a self-consistent solution compatible with the data. In this paper we examine the available RV data of the HD 181433 system (Bouchy et al. 2009) taking a more general approach, going beyond a formal fit of the Keplerian orbital elements. Even if the RV observations do not span a single period of the outermost planet, we can give reasonable constraints on the orbital elements of the poorly sampled third planet by studying the dynamics of the system. In the Keplerian fit, this important information is completely omitted.

In many situations, the interactions between planets are negligible and can be ignored. So the RV signal is just a linear superposition of different Keplerian RV curves. On the other hand, planetary interactions and resonances can be important and must be taken into account performing numerical integrations. In the short time-scale, these interactions can cause significant variations in the orbital parameters of the planets. As in our case, an ensemble of constant Keplerian orbital elements is not adequate to model the RV data and an *N*-body Newtonian model should be applied.

However, it can also happen that good Newtonian fits to the data produce planetary orbital parameters that are stable for the period of observations but lead to disruption on time-scales substantially shorter than the age of the planetary system. In this respect, long-term stability is an additional but necessary constraint that must be satisfied by multiplanet fits.

In Section 2, we perform an independent analysis of the RV data for HD 181433. We probe the phase space of the orbital parameters looking for likely configurations stable for long time-scales, say millions of years. We assume the motion is described by Newtonian interacting orbits. In Section 3, we present the dynamical study of the stable best-fitting solution and we analyse the behaviour of other plausible stable configurations. In particular, we focus on the description of secular apsidal resonances (SARs) and mean motion resonances (MMRs). In Section 4, we briefly summarize our findings, we discuss on the possibility of a terrestrial planet in the habitable zone (HZ) and we make some predictions about what we may expect from further observations.

## 2 RADIAL VELOCITY DATA ANALYSIS

Bouchy et al. (2009) announced the detection of three planets around HD 181433, a K3IV star, considering 107 RV measurements which covered more than 4 years, from 2003 June to 2008 March. The median uncertainty for the RV is 0.53 m s^{−1} with most of the data in the range 0.4–1.0 m s^{−1}. The peak-to-peak velocity variation is 48.12 m s^{−1}, while the velocity scatter around the mean RV in the measurements is 13.86 m s^{−1}.

The data do not completely cover a full period for the third planet. In fact, what is possible to spot it is just an additional long-term trend which is modelled by Bouchy et al. (2009) as being produced by a planet of minimum mass *m* sin *i* = 0.54 *m*_{Jup} on a Keplerian orbit with a period of about 6 years and *e* = 0.48. This model is unstable on the order of just thousands of years. The instability arises due to close encounters between planet *c* and planet *d*.

We perform a re-analysis of the HARPS data using the Systemic Console3 (Meschiari et al. 2009). Systemic has already been used to derive orbital fits in other works such as Vogt et al. (2010) and Meschiari et al. (2011). The list of available tools include Lomb–Scargle (LS) periodogram to identify periodicities in the RV data set, LS periodogram of residual to study periodicities in the residual RV data set, simulated annealing for global multiparameter optimization, while for local multidimensional optimization there are the Levenberg–Marquardt (LM) scheme which ensures a rapid convergence and the Nelder–Mead (sometimes called Amoeba) algorithm (Press et al. 1992).

We held the stellar mass fixed, adopting the value *M*_{*} = 0.781 M_{⊙} (Sousa et al. 2008). We make the assumption the system is coplanar and viewed edge-on. This conjecture diminishes the quantity of potential orbital configurations greatly, but the plane (*a _{d}*,

*e*) is dynamically representative for the system in the sense that it crosses all resonances (Robutel & Laskar 2001). When long enough time-series of precision data are available, the effects of mutual interactions part of the Newtonian model can potentially help in determining or estimating the masses and inclinations for the planets (see the discussion about this case in Section 2.2).

_{d}### 2.1 The Keplerian three-planet best fit

The LS periodogram shows a peak at *P* = 1171.35 d with an estimated false alarm probability (FAP) of ≈2 × 10^{−19}. Fig. 1 shows the LS periodogram of the full RV data set and the periodogram of sampling. The latter shows peaks that are related with periodicities in the cadence of observations, for instance these can arise from the solar and sidereal day, the synodic month and the solar year.

The residuals periodogram reveals an additional signal at *P* = 9.37 d with FAP of ≈3.8 × 10^{−5}. The best two-planet Keplerian fit yields residuals with an rms scatter of 2.44 m s^{−1} and reduced chi-squared χ^{2}_{red} = 15.7. The jitter for HD 181433, i.e. the jitter required to have the χ^{2}_{red} equal to 1.0, is 2.35 m s^{−1}. Fig. 2 illustrates the periodogram of residuals to the two-planet solution.

To model this long-term trend, we make the starting guess of a planet in an outer 2:1 resonance with planet *c*, we adjust the mass to match the amplitude of the signal and set a small eccentricity. At this point, a Keplerian fit using the LM algorithm naturally evolves to a solution compatible with the one by Bouchy et al. (2009). The best three-planet fit achieves a χ^{2}_{red} = 4.6 with an rms scatter of 1.34 m s^{−1} and expected jitter of 1.17 m s^{−1}. We are not aware of how Bouchy et al. have obtained a lower χ^{2}_{red} and a lower value for rms for the same solution.

The left-hand panel of Fig. 3 shows the best-fitting orbital configuration at the epoch of the first observation BJD 245 2797.87. The orbits of planets *c* and *d* cross each other, because of the occurrence of strong interactions, collisions/ejections.

### 2.2 The Newtonian three-planet stable best fit

Following the argument of Anglada-Escudé, Lopez-Morales & Chambers (2010) that eccentric orbital solutions can mimic the signal of two planets in 2:1 resonant orbits, we have also tested the hypothesis of a planet in an inner 2:1 resonance with planet *c*, but it was not possible to achieve any significant improvement to the goodness of the fit with respect to the two-planet solution.

The problem of exploring a 16-parameter phase space with stability as additional requirement can get a first simplification by arguing that the elements of the inner planets are well constrained by observations. In fact, even if we set different starting points for their parameters, the fits for them converge substantially to the same values as these signals are well sampled. A confirmation to this argument comes from the eccentricity of planet *c*, *e _{c}*, which is a very discriminating parameter towards the stability of the system: if we try to constrain

*e*to lower values the fitting we achieve is poor.

_{c}Concerning the parameters that describe planet *d*, we notice higher values for the eccentricity are preferred by the fitting. Therefore, at the end the problem can be reduced in finding for each reasonable *P _{d}* the largest value for

*e*for which planets

_{d}*c*and

*d*do not undergo instability. Likewise, we can say once we have a stable solution for a pair (

*P*–

_{d}*e*) we want to investigate if it is possible to get a different pair which generates a stable configuration having the same or lower χ

_{d}^{2}

_{red}. We can describe that as being an empirical Bayesian approach of inferring the stable best fit rather than a frequentist approach which involves a time consuming number of simulations. Here the investigation is conducted by evidences like the ones given by the collision line (see later on in this section) and the outcomes of previous simulations.

The second and third planets reside in regions spanned by a number of strong low-order MMRs (see Fig. 5). We are aware of the protective role of some MMR. For instance, the 2:1 MMR associated with the SAR consent together stable configurations even for enormously high eccentricities, ≃0.95–0.98 (see Goździewski et al. 2003, and references therein). This could explain a very large eccentricity for planet *d* and still preserve the system stability by keeping the planets away from close encounters. Actually, a modification of the relative phase of the planets strongly affects the synthetic RV curve and a *stable* resonance configuration can be far from being consistent with the RV observations. We find that manipulating the values of ω_{d} and of the mean anomaly, *M _{d}*, to get stable configurations is highly unfavourable by the RV data (i.e. poor fits are obtained). Therefore, this supports the argument arisen in the previous paragraph about performing an exploration focused on the (

*P*–

_{d}*e*) space while leaving to the algorithm the task of fitting, without constraints of any sort, the other parameters and in particular

_{d}*m*, ω

_{d}_{d}and

*M*.

_{d}To perform Newtonian orbital fits, Systemic offers different methods such as the Runge–Kutta, Hermite fourth-order and Gragg–Bulirsch–Stoer integrators. Fitting a Newtonian solution takes longer than a Keplerian model but it assures short time-scale interactions, relevant for planets *c* and *d*, are considered.

The following step is studying the stability of each distinct fit over a period of time related to the time-scale of unstable behaviours. For these long-term evolution tests, direct *N*-body integrations are applied to the orbital solutions considered. We integrate the orbits for at least 1 Myr using the Wisdom–Holman mapping integrator available in the swift software package (Levison & Duncan 1994). We use a time-step approximately equal to a twentieth (≈1/20) of the orbital period of the innermost planet. To study close encounters, we use the available Bulirsch–Stoer integrator with a tolerance parameter of 10^{−9}. We identify each configuration as being a stable system if orbits stay well bounded over an arbitrarily long period of time.

The results of our analysis are outlined in Figs 4 and 5. Here we label as stable the solutions that survived at least for 1 Myr. Once again, we fix (*P _{d}*–

*e*) and then look for the best fit, starting the LM scheme with initial points derived from previously studied configurations. The LM algorithm and Amoeba offer a clear representation of the parameter space. The dynamical analysis reveals a narrow and long band around 3.3 au and a small island around 3.2 au where good fitness is achieved and stability requirements are met. We find a configuration that we label as stable best fit, which survives for at least 250 Myr (see Section 3 for an in-depth examination). Other models scored a better χ

_{d}^{2}

_{red}but did not preserve stability for the same amount of time. Therefore, the stable best fit seems to lie on the border of a chaotic and unstable zone where small changes on the parameters of the outermost planet may push the system into a strongly chaotic state leading, in some scenarios, even to its disruption.

Fig. 4 illustrates how stable configurations do not exist in the near neighbourhood of the statistically best fit; smaller quantities for *e _{d}* are needed in order for the models to retain stability and that increases the value of χ

^{2}

_{red}.

The top panel of Fig. 5 shows the best fits obtained during our investigation in terms of the mass for planet *d*, *m _{d}*, and the semimajor axis

*a*. The picture makes it clear how to explain a certain RV amplitude

_{d}*K*, a bigger mass

_{d}*m*is required as long as

_{d}*a*increases.

_{d}The bottom panel of Fig. 5 illustrates the results of our analysis in the semimajor axis–eccentricity plane (*a _{d}*–

*e*). The parameters represented are the osculating elements at the epoch of the first observation. We show the collision line which is defined in terms of semimajor axes and eccentricities as

_{d}*a*(1 +

_{c}*e*) =

_{c}*a*(1 +

_{d}*e*). This line denotes the region where the mutual interactions of relatively massive companions can rapidly destabilize the configuration and is calculated for

_{d}*e*= 0.269 and

_{c}*a*= 1.773 au (the values are from our stable best-fitting solution; see Table 1). Note how the statistically best fit is positioned well over the collision line. We also identify the most relevant MMRs between planets

_{c}*c*and

*d*, such as the 2:1, 11:5, 9:4, 7:3, 12:5, 5:2, 8:3, 11:4, 3:1, 10:3 and 7:2. The position of the indicated locations has been calculated with respect to the values of the stable best fit. Planets in some resonant configurations, even if under the collision line, exchange angular momentum rapidly; their eccentricities are quickly pumped and that may lead again to instabilities and self-disruptions. In particular, we have found models that show a stable and bounded evolution for many Myr before the unstable behaviours are manifested. On the contrary, other resonant configurations, such as the 5:2 and 7:2, are observed to retain stability even for values over the collision line.

Parameter | HD 181433 b | HD 181433 c | HD 181433 d |

P (d) | 9.37459 | 975.41 | 2468.46 |

T_{peri} (BJD−245 0000) | 2788.9185 | 2255.6235 | 1844.4714 |

e | 0.388 40 | 0.269 12 | 0.466 26 |

ω (°) | 202.039 | 22.221 | 319.129 |

K (m s^{−1}) | 2.57 | 14.63 | 9.41 |

msin i (m_{Jup}) | 0.023 35 | 0.652 82 | 0.525 14 |

msin i (m_{⊕}) | 7.4 | 207.5 | 166.9 |

a (au) | 0.080 13 | 1.773 10 | 3.293 47 |

V (m s^{−1}) | 40 212.846 | ||

rms (m s^{−1}) | 1.36 | ||

χ^{2}_{red} | 4.96 |

Parameter | HD 181433 b | HD 181433 c | HD 181433 d |

P (d) | 9.37459 | 975.41 | 2468.46 |

T_{peri} (BJD−245 0000) | 2788.9185 | 2255.6235 | 1844.4714 |

e | 0.388 40 | 0.269 12 | 0.466 26 |

ω (°) | 202.039 | 22.221 | 319.129 |

K (m s^{−1}) | 2.57 | 14.63 | 9.41 |

msin i (m_{Jup}) | 0.023 35 | 0.652 82 | 0.525 14 |

msin i (m_{⊕}) | 7.4 | 207.5 | 166.9 |

a (au) | 0.080 13 | 1.773 10 | 3.293 47 |

V (m s^{−1}) | 40 212.846 | ||

rms (m s^{−1}) | 1.36 | ||

χ^{2}_{red} | 4.96 |

Table 1 reports the determined set of orbital elements for the stable best fit. For each planet, we list period (*P*), time of periastron passage (*T*_{peri}), eccentricity (*e*), argument of pericentre (ω), semi-amplitude (*K*), minimum mass (*m*sin *i*) and semimajor axis (*a*); we indicate also the stellar offset (*V*). This model has χ^{2}_{red} = 4.96, an rms scatter of 1.36 m s^{−1} and the expected jitter of HD 181433 is 1.19 m s^{−1}. Fig. 6 displays the RV data fitted to this model along with the residuals. The right-hand panel of Fig. 3 shows the orbital configuration of the system, this time the orbits of planets *c* and *d* do not cross each other.

Since the stable best fit is found in an active region, rather than estimating an uncertainty on each parameter, we think Figs 4 and 5 are more useful in visualizing the results of the dynamical study and highlight what is plausible to expect from new observations. If we compare our results with what has already been published for this planetary system,4 we find that the parameters of planet *b* and *c* are confirmed to be already well constrained with just *K _{c}* not compatible within the 3σ. For planet

*d*, all the elements are found within the 3σ from the original conclusion. However, it is worth to underline how to explain the very large eccentricity of the third planet and to retain a good fit to the present data, the uncertainty on the location of planet

*d*reduces dramatically to the narrow band where the 5:2 MMR is possible. Hence, this supports how a dynamical study can be fundamental in interpreting observations, producing a self-consistent model compatible with the data and giving substantial constraints on the orbital parameters.

The data do not offer any possibility of constraining the orbital inclinations. The Newtonian model cannot be particularly improved because, aside from the fact that the signal of the outer planet is not well sampled, we need to wait for secular time-scales before the variations in *i* can be spotted via the RV method. In fact, we note that for planets GJ 876 *b* and *c* which have the strongest mutual gravitational interactions, more than 11 years of observations (corresponding to more than 60 orbits of the outer planet) were used to give a reasonable estimate of the inclinations (Correia et al. 2010).

### 2.3 Additional planets?

Finally, following the argument that the proximity of the best fit to the collision line may indicate the presence of further planets (Goździewski et al. 2008), we aim to search for four-planet Newtonian solutions. The periodogram of the residuals to the three-planet solution, in Fig. 7, displays no strong peaks that would support the evidence for additional planets in the system. Apart from more distant companions, in the inner region of the system a terrestrial planet can only survive if located between planets *b* and *c*. In fact, already planet *b* is found in the proximity of the parent star and the area between planets *c* and *d* is dominated by the strong interactions that interest the two giant planets in eccentric orbits. The existence of this last planet would support the ‘packed planetary systems’ hypothesis (Barnes & Raymond 2004).

The present data do not allow making any supposition about possible outer planets (for example, an object at 7 au would have an orbital period of around 7700 d). On the other hand, with a super-Earth in the stable zone between planets *b* and *c* the fit improves. However, this signal would be at the noise level with an *F*-test of the order of 30 per cent. The *F*-test indicates the probability that a planetary model would produce a signal similar to the one due just to noise fluctuations in the data (see Marcy et al. 2005, and references therein), so additional observations are required to investigate on the presence of a super-Earth or less massive planet in this stable region.

## 3 LONG-TERM BEHAVIOUR OF THE STABLE BEST-FITTING CONFIGURATIONS

Because of the proximity of the two outermost planets, the system cannot be stable unless a resonant mechanism is present to avoid close encounters. In this section, we aim to deepen the study of the stable best-fitting configuration as well as investigate the evolution of the orbital elements, the secular resonant arguments and critical angles of some particular configurations consistent with the RV observations.

For the stable best fit, Fig. 8 shows in the subsequent panels the time evolution of the semimajor axes and of the eccentricities. Moreover, a secular resonant angle and a critical argument of the 5:2 MMR are also illustrated. Specifically, the top-left panel of Fig. 8 highlights that for the 250 Myr of the numerical integration the apocentre of planet *c* and the pericentre of planet *d* share the same region. If we get a close-in view of the situation, we notice that actually they never cross each other. In particular, the pericentre of *d* is internal to the apocentre of *c*. The former approaches the value of *a _{d}* around every 50 000 yr. It is probably a resonance that, protecting the companions from close encounters, allows the stability of the system. The top-right panel illustrates the relative large range in which

*e*and

_{c}*e*evolve. The peak-to-peak amplitude is covered in around 2500 yr only.

_{d}*e*moves in the interval 0.17–0.52 while

_{c}*e*in the range 0.17–0.50. The present eccentricities fall in the middle of these intervals, indicating that the system has been snapped in a statistically quite probable state. Also, such a large range reminds that in multiple-planet systems the orbital eccentricities can vary considerably through secular interactions on time-scales that are long compared to observational baselines but short compared to the age of the systems. Therefore, when doing statistical studies on exoplanetary systems the planetary orbits should normally be described by a complete distribution of values for the eccentricities rather than just by the present quantities (see also Adams & Laughlin 2006). This model is not observed to be in SAR. The bottom-left panel indicates the time evolution of the secular argument ω

_{d}_{c}+ω

_{d}which alternates librations with circulations. Besides, we find that 5

*n*− 2

_{d}*n*≈−3°.4 yr

_{c}^{−1}indicating the proximity of the 5

*c*:2

*d*MMR, therefore here we have a scenario similar to the Jupiter–Saturn case in the Solar system. We have studied the resonant arguments of this MMR and found that the angle 5λ

_{d}− 2λ

_{c}− 3ω

_{d}librates around 180° with a semi-amplitude of about 110°. Thus, this configuration is seen to be locked in a MMR, the bottom-right panel illustrates how this critical angle evolves with time.

The mass of planet *b* is negligible with respect to *c* and *d*, so we can assume the dynamics of the two giants are not disturbed much by the presence of the rocky planet close by. Then, we study some possible configurations with planets *c* and *d* near MMR. We do not find any plausible (χ^{2}_{red}≤ 6.03) stable solution that would correspond to the 2:1, 11:5, 9:4, 7:3, 8:3, 11:4 and 10:3 MMRs. In particular, the configurations nMMR 11:5, 9:4 and 7:3 seem preferred by the data but the excessive pumping of the eccentricities causes close encounters and planetary scatterings which do not favour stability. On the other hand, outer MMRs are possible because here in particular the planets are more spread and collisions can be avoided.

We compute the evolution of the orbital elements for fits corresponding to MMRs 12:5, 7:3, 3:1 and 7:2 which our simulations have demonstrated to preserve stability for at least 40 Myr. The results are illustrated in Figs 9, 10, 11 and 12 and show the complexity of the possible dynamical behaviours of the HD 181433 system that are consistent with the RV observations.

The case we show nMMR 12:5 has χ^{2}_{red} = 5.91 and rms scatter of 1.46 m s^{−1} with 12*n _{d}*− 5

*n*≈ 1°.8 yr

_{c}^{−1}. This model is in SAR with ω

_{c}−ω

_{d}librating around 0° with a semi-amplitude of about 45°. Moreover, we have found some critical arguments of the MMR to alternate librations with circulations implying the resonance excites a chaotic configuration. The time evolution of the critical angles 12λ

_{d}− 5λ

_{c}− 7ω

_{c}and 12λ

_{d}− 5λ

_{c}− 5ω

_{d}− 2ω

_{c}are illustrated in Fig. 9.

The scenario nMMR 7:3 has χ^{2}_{red} = 5.89 and rms scatter of 1.46 m s^{−1} with 7*n _{d}*− 3

*n*≈ 8°.7 yr

_{c}^{−1}. This model is observed to be in SAR with the critical angle ω

_{c}−ω

_{d}librating around 0° with a semi-amplitude of about 40°, meaning that their periastrons are aligned. We have studied the resonant arguments of the 7:3 MMR and found that three angles, i.e. 7λ

_{d}− 3λ

_{c}− 4ω

_{d}, 7λ

_{d}− 3λ

_{c}− 3ω

_{d}−ω

_{c}and 7λ

_{d}− 3λ

_{c}− 2ω

_{d}− 2ω

_{c}, alternate librations with circulations. This indicates that the configuration is close to the resonance separatrices. Fig. 10 illustrates how the critical angle 7λ

_{d}− 3λ

_{c}− 3ω

_{d}−ω

_{c}evolves with time.

The case near the low-order MMR 3:1 has χ^{2}_{red} = 6.03 and rms scatter of 1.47 m s^{−1} with 3*n _{d}*−

*n*≈ 0°.6 yr

_{c}^{−1}. This model is in SAR with ω

_{c}−ω

_{d}librating around 180° with a semi-amplitude of about 110° implying that this time the periastrons of planets

*c*and

*d*are anti-aligned. The data points that diverge from the periodic signal (Fig. 11, top-right panel) represent the instants when

*e*gets close to be null and so the argument of pericentre, ω

_{d}_{d}, is not well defined. Furthermore, for the critical arguments of the MMR, we find that librations alternate with circulations indicating a chaotic zone spanned by overlapping resonances. The time evolution of the critical angles 3λ

_{d}−λ

_{c}−ω

_{d}−ω

_{c}and 3λ

_{d}−λ

_{c}− 2ω

_{c}are also illustrated in Fig. 11.

The scenario nMMR 7:2, illustrated in Fig. 12, has χ^{2}_{red} = 5.65 and rms scatter of 1.45 m s^{−1} with 7*n _{d}*− 2

*n*≈−1°.7 yr

_{c}^{−1}. This model is not seen to be in SAR. However, we find it to be locked in the MMR as the critical argument 7λ

_{d}− 2λ

_{c}− 5ω

_{d}librates around 180° with a semi-amplitude of about 85°. It is worth notice how this mechanism is capable of pumping

*e*from 0.09 to 0.47 in less than 5000 years. This configuration is located over the collision line; it is the resonance that prevents close encounters and provides long-term orbital stability to the system.

_{c}Considering all the configurations studied, the behaviour of planet *b* seems to be unrelated to the two giant companions as the amplitude of the eccentricity signal of *b* appears to be unaffected even in the cases in which *c* and *d* are trapped in a resonance. Moreover, it seems unrealistic that planet *b* can be involved in a *p*:*q*:*r* MMR with the outer two planets: *b* is too far away from them and, in addition, it is mainly influenced by general relativistic and tidal effects (a discussion of these mechanisms goes beyond the aim of this paper).

The synthetic RV curves for the Keplerian best fit, stable best fit and the models nMMRs 12:5, 7:3 and 7:2 are illustrated in Fig. 13. The period through the year 2017 is covered. It is difficult to distinguish the curves from each other in the time range covered by the observations, but at the time of writing the difference can be spotted. However, the signal of the Keplerian best fit will diverge evidently from the stable best fit only in 2013 February (∼JD 245 6340).

## 4 DISCUSSION AND CONCLUSION

In this paper, our efforts have been focused on finding a plausible solution to the available RV data for the planetary system of HD 181433. In our investigation, we have included an analysis of the long-term evolution of the system, and the results support the thesis that the dynamics is an important observable that has to be taken into account with the same priority as the RV observations.

The story with HD 181433 does not differ from the one of many other multiplanet systems found on the edge of long-term dynamical stability. The dynamical modelling of the RV with stability constraints offers precious information on the dynamical architecture of the putative planetary configurations. The stability criterion becomes a fundamental tool which provides limits to the physical and orbital elements of the planets when data do not cover completely the longest orbital period.

Our investigation leads to a Newtonian model for HD 181433 stable for at least 250 Myr. The solution is compatible with what was found by Bouchy et al. (2009), but our analysis strongly diminish the uncertainty on the location of planet *d* to the exiguous band where the 5:2 MMR is possible and stability is preserved. This seems the only plausible way to explain a very large eccentricity for the outermost planet, a quality which must be met in order to hold a good fit to the present data. In general, we can say that when an unstable high eccentric solution is found for a multiplanet system, the study of resonances may lead to the finding of a reasonable stable solution. By doing so, it is possible to constrain with high confidence the orbital period of the outermost (poorly sampled) planet well before sufficient data, covering several orbital periods, become available.

Apart from the 5:2 MMR, the orbital evolution of the two giant companions is confined to a zone spanned by a number of other low-order two-body MMRs. We have studied different plausible stable configurations for the planetary system and, in particular, we have illustrated the behaviours caused by SARs and MMRs. We have also found that at the time of writing with new data points it is definitely feasible to refine the circle of likely scenarios.

Furthermore, given the strong gravitational interactions between the two giant planets, a self-consistent *N*-body model for the RV data will help in estimating the inclination of the planetary orbits and of the physical masses (the RV method returns just the minimum masses for the planets). If it is not possible to use in synergy other methods, e.g. astrometry, transits, it may be necessary to observe around 50 full orbits of planet *d*, i.e. 300 years of RV data, to strongly constrain the orbital inclinations. As the large values already observed for the eccentricities may have been trigged by mechanisms which influence also the orbital inclinations (e.g. Libert & Tsiganis 2009), a considerable value for their mutual inclination may be expected.

We can calculate the HZ orbital distance *a*_{hab}, defined as the distance where a planet would receive the same insolation as the Earth: au. For HD 181433, we find *a*_{hab} = 0.55 au (*P*_{hab}∼ 170 d) which is in the region between planets *b* and *c*. Our simulations, run for 40 Myr, shows that an Earth-size planet in the HZ (and in eccentric orbit) can retain stability indeed. The existence of a planet in this zone not only fills an empty gap in the system, but would also offer a harbour for life. Therefore, this hypothetical planet becomes interesting for a double reason. Additional observations are required to investigate on the presence of further bodies and in particular on the existence of a terrestrial planet in the HZ of HD 181433.

The planetary system of HD 181433 offers a wide range of challenges that can go from understanding its real physical architecture to the study of potentially habitable worlds. Further observations can confirm the results illustrated in this paper, improve our understanding of the dynamical structure of this system and, in general, give additional insights into the study of the dynamics of planetary systems.

The author acknowledges support of studentship from Queen Mary University of London and would like to thank Craig Agnor, Carl Murray, Stefano Meschiari and Richard Nelson for useful discussions and suggestions while conducting this research. We are also grateful to the anonymous referee for the valuable comments which improved the results presented in this paper.