© ESO, 2015

1. Introduction

With radial velocity (RV) measurements of the stellar reflex motion caused by orbiting companions a pair of planets in low-eccentricity orbits near a 2:1 mean motion resonance (MMR) can be misinterpreted as a single planet with moderate eccentricity (Anglada-Escudé et al. 2010; Wittenmyer et al. 2013). This is, in particular, possible when the available data are sparse and have large errors, when the stellar RV amplitude induced by the inner planet is smaller than that induced by the outer planet, and when the overall scatter of the data around either model is too large to make it possible to distinguish them.

A one-planet system might intuitively be considered a simpler model than a two-planet system and would therefore be favoured when applying Occam’s razor to select the simplest hypothesis. However, pairs of planets near a 2:1 MMR are not rare and so warrant consideration. Lissauer et al. (2011) find that at least 16% of the systems with more than one candidate for a transiting planet in the data from the Kepler satellite mission (e.g. Borucki et al. 2010) include a pair of planets with period ratios in the range 1.82:1−2.18:1, i.e. within 9% of the “pure” 2:1 ratio. Delisle & Laskar (2014) find an excess of planets with a period ratio a few percentage points higher than the 2:1 or 3:2 resonant value in data from the Q1-Q16 KOI catalogue (Batalha et al. 2013). This effect is most pronounced for systems with periods of the inner planet < 5 d, for which no systems in the exact resonance are actually found, it is still significant for inner planet periods between 5 d and 15 d, but strongly reduced for longer periods, for which the exact resonance is observed more frequently. In their more recent analysis of Kepler data, Steffen & Hwang (2015) find a significant tendency for pairs of planets to have a period ratio close to 2.2:1. In RV data, the true period ratio will be difficult to distinguish from the MMR when data are sampled sparsely and when they only cover a few planetary orbits.

Furthermore, it is interesting to note that the model of a two-planet system in precisely 2:1 resonant circular orbits (hereafter abbreviated as “2:1-RCO”) has the same number of free parameters as the model of a single planet in a Keplerian eccentric orbit. In the case of the 2:1-RCO, the stellar RV varies as a function of time t as $\begin{array}{ccc}\mathrm{RV}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \mathit{\gamma }\mathrm{+}{\mathit{\kappa }}_{\mathrm{1}}\cos \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\max \mathit{,}\mathrm{1}}}{\mathit{P}}\right)\mathrm{+}{\mathit{\kappa }}_{\mathrm{2}}\cos \left(\mathrm{4}\mathit{\pi }\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\max \mathit{,}\mathrm{1}}}{\mathit{P}}\mathrm{-}\mathrm{\Delta }\mathit{\phi }\right)\\ & \mathrm{=}& \mathit{\gamma }\mathrm{+}{\mathit{\kappa }}_{\mathrm{1}}\cos \left(\mathrm{2}\mathit{\pi }\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\max \mathit{,}\mathrm{1}}}{\mathit{P}}\right)\mathrm{+}{\mathit{\kappa }}_{\mathrm{2}}\cos \left(\mathrm{4}\mathit{\pi }\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\max \mathit{,}\mathrm{2}}}{\mathit{P}}\right)\mathrm{·}\end{array}$Here P is the period of the outer planet, κ₁ and κ₂ are the RV semi-amplitudes of the star induced by the outer and inner planet, respectively, γ is the systemic RV, and Δφ is the phase difference between the two sinusoids. Finally, t_max,1 and t_max,2 are the times of the RV maximum that would occur if there were only the outer, respectively, the inner planet in the system¹.

In the case of a single planet in an eccentric orbit the stellar reflex motion is given by $\mathrm{RV}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{\gamma }\mathrm{+}\mathit{K}\frac{\sin \mathit{E}\sin \mathit{\omega }\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\cos \mathit{E}\cos \mathit{\omega }}{\mathrm{1}\mathrm{-}\mathit{e}\cos \mathit{E}}\mathrm{·}$(3)Here E is called the “eccentric anomaly”, which is related to the “mean anomaly” M via the Kepler equation $\mathit{E}\mathrm{-}\mathit{e}\sin \mathit{E}\mathrm{=}\mathit{M}\mathrm{=}\mathrm{2}\mathit{\pi }\frac{\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{p}}\mathrm{)}\mathrm{mod}\mathit{P}}{\mathit{P}}\mathrm{·}$(4)After introducing the “true anomaly” θ given by $\tan \mathit{\theta }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\frac{\sin \mathit{E}}{\cos \mathrm{(}\mathit{E}\mathrm{-}\mathit{e}\mathrm{)}}\mathit{,}$(5)the expression for RV(t) can be written as $\mathrm{RV}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{\gamma }\mathrm{+}\mathit{K}\mathrm{\left[}\cos \mathrm{\right(}\mathit{\theta }\mathrm{+}\mathit{\omega }\mathrm{)}\mathrm{+}\mathit{e}\cos \mathit{\omega }\mathrm{]}\mathit{.}$(6)The free parameters of this Keplerian model are the orbital period P, the stellar RV semi-amplitude K, the orbital eccentricity e, the longitude of periastron ω, the time of periastron passage t_p, and the systemic RV γ.

Both models have six free parameters. It is a known difficulty with Occam’s razor that the term “simplicity” is not easy to define. On the one hand, the mathematics behind the Keplerian orbit is more complicated than that of circular orbits, since one has to solve the transcendent equation for the eccentric anomaly E (Eq. (4)); also the description of the Keplerian orbit (Eqs. (3) and (6)) involves more trigonometric functions than that of the 2:1-RCO pair of planets (Eqs. (1) and (2)). On the other hand, the number of assumptions leading to Eqs. (1) or (2) is larger than for the single eccentric planet, since the eccentricities of both planets of the pair must (for practical purposes, i.e. within the measurement errors) be exactly zero and the period ratio must be exactly 2:1. If these assumptions are not made, and in reality they cannot be fulfilled exactly, the pair of planets has more model parameters than the single planet model, namely eleven vs. six for the case that the system is well described by two Keplerians. In general, however, a realistic model for a pair of planets is even more complicated when the gravitational interaction between the two planets is considered. Their mutual gravitational disturbance leads to evolving orbital parameters (e.g. Wittenmyer et al. 2013 and references therein; Trifonov et al. 2014) which are unlikely to be found occupying the exact values required for the 2:1-RCO case at any time.

The discussion which model is simpler and therefore to be preferred has, of course, a precursor in history important for our overall worldview, namely the question whether Kepler’s elliptical orbits were to be preferred over Ptolemy’s epicycles when it came to describing the observed positions of those of the solar system planets known at the times of these historical astronomers².

Resonant orbits are found in the solar system as well as in extrasolar planetary systems. The most famous solar system examples for 2:1 resonant orbits are Jupiter’s moons Ganymede, Europa, and Io which are in a configuration called the Laplace resonance with orbital period ratios very close to 4:2:1. Other examples for 2:1 resonances in the solar system are Saturn’s moons Tethys and Mimas as well as Dione and Enceladus. Also, the orbital period ratio of Neptune and Uranus deviates by less than 2% from the 2:1 resonance. As an example for the possible destabilizing effect of resonances the Cassini division in Saturn’s rings should be mentioned which has been cleared by a 1:2 resonance with the moon Mimas. Concerning the fact that resonances with massive bodies can stabilize or destabilize orbits we note that, on the one hand, the Trans-Neptunian Objects display various types of resonances with Neptune among which also the 2:1 ratio is found in several cases; on the other hand 1:2 resonances of Main-Belt asteriods with Jupiter are very rare.

The Laplace resonance is also found for the components e, b, and c of the planetary system around the M4V dwarf GJ 876 (Rivera et al. 2010; Marcy et al. 2001) consisting of two Jovian planets in a 2:1 resonance (the more massive one being in the outer orbit of the pair) plus a third low-mass planet in the 4:2:1 resonant orbit further out. The system has a fourth, non-resonant Super-Earth-type planet in a very tight inner orbit around the host star. Further examples for 2:1 resonant pairs of planets include the system around the G0V star HD 82943 with two Jovian planets (Mayor et al. 2004) of similar minimum mass. Indications presented by Beaugé et al. (2008) for a third lower-mass planet in an outer orbit and in a 4:2:1 resonance with the inner two were refuted by Tan et al. (2013). A pair of two Jovian planets near the 2:1 resonance was also found in the system of the K0V star HD 128311, the outer one being more massive. There is also evidence for a third inner non-resonant Saturn-type planet (McArthur et al. 2014; Vogt et al. 2005; Butler et al. 2003). Remarkable in this system are the relatively high eccentricities of 0.303 and 0.159 for the inner and outer planet, respectively, of the resonant pair. Around the K-giant η Cet, Trifonov et al. (2014) found a system consisting of two Jovian planets whose period ratio is within ≈ 9% of the 2:1 resonance. Another special system is KOI-730 whose four supposed Super-Earth planets seem to have period ratios of 3:4:6:8, which thus include two interlaced occurences of the 2:1 resonance (Lissauer et al. 2011).

In the present paper we investigate the possibility that some fraction of the extrasolar planets in eccentric orbits found by the RV technique are actually pairs of planets in orbits that are close to circular and near a 2:1 resonance. By simulating and comparing the respective orbits as a function of their model parameters (Sects. 2 and 3) we explore the circumstances under which the two models are indistiguishable. This enables us to identify those of the known eccentric planets in the literature for which suitable follow-up observations could determine which model is correct. To illustrate this, we study one such example (Sect. 4). In the following discussion we provide recommendations for the strategy for such new observations that may uncover a so far unknown additional planet near the 2:1 resonant orbit (Sect. 5). Finally, we summarize our main conclusions (Sect. 6).

2. Simulations

We have simulated RV data for pairs of 2:1-RCO planets for the purpose of fitting them with single-component Keplerian models as well as with single-component circular models (i.e. sine waves). For this we used the fitting routines provided by the software for the generalized Lomb Scargle (GLS) periodogram (Zechmeister & Kürster 2009). The simulated data RV_i are composed of two components RV_1,i and RV_2,i and calculated from $\begin{array}{ccc}{\mathrm{RV}}_{\mathit{i}}\mathrm{=}{\mathrm{RV}}_{\mathrm{1}\mathit{,i}}\mathrm{+}{\mathrm{RV}}_{\mathrm{2}\mathit{,i}}\mathrm{=}{\mathit{\kappa }}_{\mathrm{1}}\cos \mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}{\mathit{\kappa }}_{\mathrm{2}}\cos \mathrm{(}\mathrm{4}\mathit{\pi }{\mathit{t}}_{\mathit{i}}\mathrm{-}\mathrm{\Delta }\mathit{\phi }\mathrm{)}\mathit{.}& & \end{array}$(7)A comparison with Eq. (1) shows that we have chosen the period of the outer planet as P = 1 (hence the period of the inner planet is 0.5) and set the uninteresting parameters to t_max,1 = 0, and γ = 0.

2.1. Dense sampling

As a numerical approximation to continuous RV curves, we created densely sampled data sets containing 1000 points calculated at equidistant sampling times t_i and filling the interval [ 0,P), i.e. the sampling times are given by the values 0.000, 0.001,...,0.999 and provide an even distribution of the sampled phases. We produced these data with the following 17 values of the ratio of semi-amplitudes: κ₂/κ₁ = i/8 for i = 0,...,16, i.e κ₂/κ₁ = 0, 1/8, 1/4, 3/8, ..., 2. The semi-amplitudes in Eq. (7) were always normalized according to $\mathit{\kappa ̂}\mathrm{1}\mathrm{=}{\mathit{\kappa }}_{\mathrm{1}}\mathit{/}\mathit{\kappa }\mathrm{and}\mathit{\kappa ̂}\mathrm{2}\mathrm{=}{\mathit{\kappa }}_{\mathrm{2}}\mathit{/}\mathit{\kappa }\mathrm{with}\mathit{\kappa }\mathrm{=}\sqrt{{\mathit{\kappa }}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\kappa }}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}$(8)where and denote normalized values, whereas κ₁ and κ₂ stand for non-normalized values.

The simulated data were subsequently fitted with four different single-component models:

• 1.

  sine fit with a period of P = 1,

• 2.

  sine fit with a period of P = 0.5,

• 3.

  Keplerian fit with a period of P = 1,

• 4.

  Keplerian fit with a period of P = 0.5.

For each value of the semi-amplitude ratio of the inner over the outer planet we produced 12 data sets with the 12 different, equally spaced phase shift values Δφ = 0°, 30°,...,330°. A natural exception to this is the data set with for which the phase shift is not defined.

In order to compare the fits, their rms deviations σ_model from the simulated 2:1-RCO data were employed to indicate the quality of the fit. This is justified since our simulations are based on idealized data without errors, and we compare the optimum fits of either type without any uncertainties in the fit parameters. However, whenever we draw any conclusions from real data (Sects. 4 and 5) we instead compare the probability of chi-square p(χ²) of the fits. We express rms values for both the Keplerian fits and the sine fits in units of the semi-amplitude K obtained for the Keplerian fit with P = 1.

Fig. 1 Examples for single-component fits to the RV data of a pair of 2:1-RCO planets. Thick solid lines represent the input data which are shown for amplitude ratios and 2/1 (from left to right) and values of the phase shift Δφ = 0° (upper row of panels) and 90° (lower row of panels). Thin solid lines depict the best-fitting Keplerian orbits with P = 1 (not discernible in the leftmost panels because of their very good match to the thick solid line). Dotted lines indicate best-fitting circular orbits (sine waves) with P = 1. Dashed lines stand for best-fitting circular orbits with P = 0.5. Note that the best-fitting Keplerian orbits with P = 0.5 are identical to the circular orbits with the same period, since the eccentricity vanishes for all fits. Plots for phase shifts of 180° and 270° can be obtained by flipping the 0° and the 90° plots, respectively, around the horizontal axis and then shifting them by 0.5 in phase.

Fig. 2 Results from our simulations of fitting 2:1-RCO data with various models. Left panel: the rms deviation σmodel between the various types of fit and the 2:1-RCO input data as a function of the amplitude ratio of the input data. Shown is σmodel for the single-component Keplerian fits (“KF”) and for the sine fits (“SF”) pertinent to the circular models. Note that σmodel is expressed in units of the semi-amplitude K of the Keplerian fit with Period P = 1. Middle panel: variation of the eccentricity e of the Keplerian fits with . Right panel: the semi-amplitude K of the Keplerian fits as a function of and normalized to one for . In all panels fits with period P = 1 are plotted with larger symbols than fits with P = 0.5 (see the labels at the curves). Thin dotted lines in all three panels: as a comparison we have included the result from a simulation for sparse data sampling based on the times at which the 20 RV measurements for the host star of HD 27894 were taken by Moutou et al. (2005). These thin dotted lines correspond to the median values for σmodel, e, and K, respectively, obtained from this simulation (Keplerian fits with P = 1 only). Note that for less than ≈ 0.5 they are located near the line for KF (P = 1) for e and K, whereas for σmodel the similarity holds up to .

2.2. Sparse sampling

In order to compare the results from our idealized densely sampled data sets with those that are obtained from real observations, we created rather sparsely sampled data sets consisting of only 20 data points sampled at the temporal pattern of the observations of the planet host star HD 27894 (Moutou et al. 2005). These observations are quite unusual due to the small number of data points that led to the discovery of this planet.

With this sparse sampling grid we simulated again 2:1-RCO RV data with the same 17 values for the amplitude ratio κ₂/κ₁ and the same 12 phase shift values Δφ as for the densely sampled data set described in Sect. 2.1. We then fitted these data, but restricted ourselves for the purpose of this comparison to only the Keplerian model with period P = 1. Due to the larger number of free parameters compared with the sinusoidal fits the strongest discrepancies are expected for the Keplerian fits.

3. Results

3.1. Single-planet fits to 2:1-RCO planets

Figure 1 shows examples of single-planet model fits to the densely sampled simulated input data with different amplitude ratios and two selected values of the phase shift Δφ = 0° and 90°, while Table 1 shows the mean rms deviation σ_model of these model fits from the data where the mean has been taken over all 12 phase shift values (see Sect. 3.2 on the relatively small variation of the individual values of σ_model as a function of the phase). Actually, as can be seen from Fig. 2 (left panel), the single Keplerian always yields the best fit, i.e. the smallest σ_model, for , while for the fit with P = 0.5 is better.

Specifically, for an amplitude ratio (Fig. 1, left panels), the data (thick solid line) and the Keplerian fit with the period of the outer planet (thin solid line) become indistiguishable in our plot, whereas the single sinusoidal fit with the same period (dotted line) can be seen to systematically deviate from the data, and the fit with the period of the inner planet (dashed line) is completely inadequate. For an amplitude ratio of 1/2 (Fig. 1, second column of panels), the single Keplerian with the period of the outer planet deviates systematically from the data, but could still be a reasonable match to real data, if signal-to-noise ratios are moderate. For this amplitude ratio, the sinusoidal fit with the period of the outer planet deviates much more from the data, and the fit with the period of the inner planet deviates even more (but less than for the amplitude ratio of 1/8). For an amplitude ratio of 1/1 (Fig. 1, third column of panels), all single-planet fits are inadequate, but the Keplerian fit with the period of the outer planet is still the formally best fit. For this amplitude ratio, the sine fits with P = 1 and P = 0.5 are of equally poor quality since they have the same rms deviation (see Table 1). If the amplitude ratio of the inner vs. the outer planet is reversed then for the fit with the period of the inner planet is the best one, followed by the Keplerian fit with the period of the outer planet, and then the sine fit with the period of the outer planet (see Fig. 1, right panels which are for ).

3.2. Limited dependence on phase shift

For the densely sampled data sets, we find that the interesting fit parameters RV semi-amplitude K and eccentricity e, and also the rms of the fit vary as a function of phase shift Δφ only by very small amounts³. For the extreme amplitude ratio , these variations are 0.55%, 1.5%, and 8.0% peak-to-valley, respectively, for the rms, for K, and for e. For smaller amplitude ratios, these values are even smaller. Therefore, in our analysis, for every given amplitude ratio, we are using the mean values for the rms and for e and K obtained from the fits to the 12 data sets with different phase shifts (see Sect. 3.1). For the sparsely sampled data sets, we use the median values (instead of the mean values) of e and K for all further comparisons, since they display somewhat stronger variations and can have strong outliers so that the median is the more representative value.

3.3. Zero eccentricity in P = 0.5 fits

We also find that in the case of the Keplerian models with period P = 0.5 all fits to the densely sampled data sets yield an eccentricity e = 0 and an amplitude and rms equal to that of the sine fits with the same period value. That is why we will henceforth no longer distinguish between Keplerian fits and sine fits for P = 0.5. As noted above, we have not made any fits with P = 0.5 to the sparsely sampled data sets. In this case, the Keplerian fit with P = 0.5 will normally show differences from the circular fit, since the sparse sampling does no longer guarantee the symmetric distribution of the RVs which is naturally present for a densely sampled sine wave. However in this paper, we will largely concentrate on cases of small deviations from the single sinusoidal model, such as low-eccentricity Keplerians or small-amplitude second circular planets with half the period of the (circular) first one. Consequently, we are exploring the regime where P = 0.5 fits alone are never a good match anyway.

3.4. Dependence of rms, e, and K on amplitude ratio for dense sampling

This behaviour can also be deduced from Fig. 2 which shows the dependence of the rms of the different types of fits on the amplitude ratio of the input data (left panel). Figure 2 also shows the variation of the eccentricity e (middle panel) and the semi-amplitude K of the Keplerian fit (right panel) with this ratio.

3.5. Example for dependence of rms, e, and K on amplitude ratio for sparse sampling

In the three panels of Fig. 2, thin dotted lines show the results from the data sets simulated using the sparse sampling example outlined in Sect. 2.2. With this example, we illustrate that conclusions concerning the goodness-of-fit, the determined eccentricity, and the RV semi-amplitude drawn from our idealized densely sampled data are likely to hold also for realistic data sets as long as the RV amplitude ratio of the 2:1-RCO planets is not too large. In particular, our sparse data example and our dense data sets produce similar σ_model values for amplitude ratios up to ≈ 1.4, and similar e and K values for up to ≈ 0.5 (see Fig. 2).

Fig. 3 Rms deviation σmodel of the Keplerian fit with P = 1 from the 2:1 resonant input model (solid line with points and labels). This deviation is expressed in units of the RV semi-amplitude K and plotted as a function of the eccentricity e determined by the fit. The labels indicate the amplitude ratio of the input sinusoids. Thin dashed line: amplitude of the first harmonic in the Fourier expansion of the Kepler equation.

3.6. Dependence of rms on e – a first tool to search for erroneous models

Figure 3 shows the rms deviation of the Keplerian fit with P = 1 from the simulated 2:1-RCO data as a function of the pertinent value for the eccentricity obtained by the same fit. This is based on the densely sampled data. We will use this relation as a tool to identify those of the known RV determined planetary systems for which single objects on eccentric orbits have been fitted, but a pair of 2:1 resonant planets is also a possible model.

Looking, in observed data (with not too sparse sampling and moderate amplitude ratios; see Sect. 3.5), at the scatter σ_Kep of the RV measurements around the fitted Keplerian model (expressed in units of the RV semi-amplitude K of the fit), Fig. 3 allows us to determine the expected (typical, average) contribution σ_model to the scatter arising from the wrong model fit. This means that all systems that we can identify in the literature which have a planet in a putative eccentric orbit as well as a scatter of the data around the eccentric model larger than the one given by the solid curve in Fig. 3 are potential candidates for a 2:1 resonant pair of planets.

In principle, the instrumental noise σ_inst and the RV jitter σ_star arising from stellar effects such as activity, convective motions, and pulsations, commonly named “stellar jitter”, should also be added in quadrature to the scatter due to the wrong model fit. Systems with an observed RV scatter σ_Kep exceeding this combined scatter would then be candidates for a follow-up study in search for a 1:2 resonant second planet, i.e. sytems with ${\mathit{\sigma }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\ge }{\mathit{\sigma }}_{\mathrm{model}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{inst}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{star}}^{\mathrm{2}}$. In practice, however, this can be done only for those systems for which good estimates of the instrumental noise and of the stellar noise are available. Estimates of the instrumental noise can often be obtained from a joint analysis of the data of the ensemble of stars observed in a survey-type programme, e.g. by looking for the most constant stars, but even for those substantial effort is required to disentangle the instrumental noise from stellar jitter. Short-term jitter due to stellar p-mode oscillations can be estimated using the scaling relation from Kjeldsen & Bedding (1995). There is also quite some literature on estimates of the effects of stellar activity (e.g. Boisse et al. 2011; Saar & Donahue 1997; Saar et al. 1998; Saar & Fischer 2000), but because of the nature of this activity (spot migration, activity cycles, etc.) they cannot be very precise.

3.7. Probability of chi-square – the tool of choice for real data

As already alluded to in Sect. 2.1, it should be noted that the scatter alone is not a statistically sound criterion, since in general, it is possible to find signals that are smaller than the scatter/noise, if one has a sufficiently large number of measurements. Rather than looking at the rms, the appropriate approach would be to look at the probability of χ² of the fitted models, and in fact our main data base for planetary orbital parameters (the EOD; see Sect. 4) provides values of χ² (along with the also required values for the degrees of freedom, d.o.f.) for a good fraction of the interesting entries, but only so for the Keplerian fit, not for an alternative fit with a 2:1 resonant pair of sinusoids.

3.8. First harmonic of the Fourier expansion of the Kepler equation

Anglada-Escudé et al. (2010) base much of their study of 2:1 resonant vs. Keplerian models on the Fourier expansion of the Kepler equation. While we have adopted a different approach in our simulations, a comparison is in order at this point. It is given by the dashed line in Fig. 3 which represents the amplitude of the first harmonic of the Fourier expansion of the Kepler equation (from Anglada-Escudé et al. 2010 and Lucy 2005). The similarity between the two curves in Fig. 3 demonstrates that the ability to distinguish between the Keplerian model and the model with the 2:1 resonant pair of circular orbits is largely equivalent to the ability to detect this first harmonic.

3.9. The model-related scatter in real data

We also note that the so determined model-related contributions σ_model to the scatter are exact only for very densely sampled RV curves such as those we have used for our simulations (Sect. 2.1). For much more sparsely sampled real data, the model-related contributions to the scatter can deviate somewhat in the individual case from the here determined values. Typically, data sets with fewer measurements can be fitted “better”⁴ by any multi-parameter model, so that the model-related contribution to the rms scatter is likely to be somewhat smaller than what we have determined for our simulated data sets with 1000 equidistant data points. Nevertheless, we will conservatively use these larger values to guide us in our search for candidate systems (i.e. systems with even larger scatter).

4. Follow-up candidates from the literature

4.1. Sample definition

In order to come up with a compilation of candidate systems, we used the Exoplanet Orbit Database (EOD)⁵ described by Wright et al. (2011). From the list of detected exoplanets, we have selected two versions of a sample of stars based on the following criteria.

• 1.

  Stars with an RV orbital solution.

• 2.

  Stars for which no planetary transits have been observed.

• 3.

  e> 0.

• 4.

  Listed number of components equal to one (supposedly single planets).

• 5.

  K must be given.

• 6.

  a) Values for the rms scatter of the data around theKeplerian fit σ_Kep are available. b) Values for the probability of chi-square $\mathit{p}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\right)}$ can beobtained from the data base as it provides both the re-duced χ² and the number of observations.

As of January 13th, 2015, we find that 254 of the stars contained in the EOD fulfill criteria 1−5 plus 6a; these will make up our version 1 sample. 187 of the stars fulfill criteria 1−5 plus 6b and will constitute our version 2 sample. 185 stars fulfill criteria 1−5 plus 6a and 6b. Concerning criterion 6b, we note that the EOD contains entries for the (reduced) chi-square ${\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}$ and the number of measurements N from which we determine $\mathit{p}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\right)}$ using standard routines (e.g. Bevington 1969)⁶.

Note that in order to constrain the sample, we discarded known transiting systems (criterion 2), since for those chances are higher that a second planet further inside in the system has already been detected through its own transits. We also restricted ourselves to supposedly single-planet systems (criterion 4), because finding an additional low-amplitude planet in a 2:1 resonance with one of the planets in the system will be more difficult when more than one planet has already been found to contribute to the observed stellar RV signal.

Fig. 4 Rms scatter σKep of the RV solutions around single-planet Keplerian fits for our sample of 254 stars for which the EOD provides these values (dots). The scatter values are plotted as a function of the formal eccentricity of the Keplerian fit and in units of its RV semi-amplitude. 1σ-errors of the eccentricity are represented by grey horizontal bars. In two instances no errors are available; the pertinent data points are marked by small circles around them. The large circle and cross marks the case of HD 27894 that is studied in Sect. 4.3. The solid curve is identical to the solid line in Fig. 3 representing the deviation σmodel between the single-planet Keplerian model and the model for the 2:1 resonant pair of planets. Scatter values above the curve are larger than the expected deviation arising from applying the erroneous model (single Keplerian instead of 2:1 resonant pair) so that both models are possible, whereas scatter values below the curve are too small to assume an erroneous model and therefore favour the single Keplerian.

Fig. 5 Cumulative histogram of $\log \mathrm{\left[}\mathit{p}\mathrm{\right(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}$ of the Keplerian fits for those 187 stars in our sample for which both the reduced chi-square ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ of the fit and the number of measurements is provided by the EOD. For a substantial fraction of the stars, the chi-square probability is too small to accept the single Keplerian model, e.g. 39% of the stars have $\log \mathrm{\left[}\mathit{p}\mathrm{\right(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\mathit{<}\mathrm{-}\mathrm{3}$ (see also Table 2).

Fig. 6 Fit to the HD 27894 RVs with a model consisting of two planets in circular orbits. Note that in this case the exact 2:1 resonance is not forced, but the two period values are both free parameters of the fit. Left: time series together with the best-fitting two-planet circular model (top, thick solid line) and residuals (bottom). For display purposes the x-axis has been broken between BJD – 2 452 950 and BJD – 2 453 267.5 where there is a large data gap. The thin dashed line in the upper panel represents the best-fitting single Keplerian model for comparison. Right: RV data phase folded with the period of the outer planet after subtraction of the signal of the inner planet (top) and phase folded with the period of the inner planet after subtraction of the signal of the outer planet (bottom).

4.2. Resulting candidate sample

Figure 4 shows the location in the rms-vs.-eccentricity diagram of our version 1 sample stars in comparison to the curve from Fig. 3, i.e. the curve that depicts the expected contribution to the rms if a 2:1 resonant system with two circular orbits is fitted with a single-star Keplerian model. We find that 188 of the 254 stars in the sample are located above the curve, i.e. 74% of the stars are possible candidates for hosting a system of two planets in a 2:1 resonance when the criterion of exceeding the model-related scatter for 2:1-RCO systems is applied (see Sect. 3.6).

In order to apply also the probability of χ² of the single Keplerian fits as a criterion for identifying candidates (Sect. 3.7), we display in Fig. 5 (on a logarithmic abscissa) the cumulative histogram of $\mathit{p}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\right)}$ for our version 2 sample stars which demonstrates that a large fraction of the stars has a very low probability that the applied model alone is correct. This is also shown by Table 2 that lists how many stars have $\mathit{p}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\right)}\mathit{<}\mathrm{0.05}\mathit{,}\mathit{<}\mathrm{0.01}\mathit{,}\mathit{<}\mathrm{0.005}$, and < 0.001, respectively, and which fraction of the total of 187 stars in our version 2 sample this corresponds to. It also provides the confidence for the rejection of the single Keplerian model for these stars given by $\mathrm{1}\mathrm{-}\mathit{p}\mathrm{(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}$). The hypothesis that the Keplerian model alone can explain the observed variability can be rejected for 54% of our version 2 sample stars at the 95% confidence level and for 39% at the 99.9% confidence level.

4.3. A case study: HD 27894

From our (version 1) sample of candidate 2:1-RCO systems (systems above the curve in Fig. 4), we selected the 20 very sparsely sampled RV measurements of HD 27894 (Table A.2 in Moutou et al. 2005) in order to provide an example for a detailed study of a supposedly single-planet system which might hide a second planet in an inner orbit near a 2:1 MMR with the first one. (We have already employed the observing times of this star in Sect. 2.2.) HD 27894 is a low-activity K2V star with a mass of 0.8 M_⊙. Further characteristics of the star are described in Moutou et al. (2005).

HD 27894 was selected because the moderate eccentricity of the published Keplerian fit as well as the considerable scatter of the data around it made this system a perfect candidate for an ambiguous, undecided case well suited for an illustrating example. This system earns much of this ambiguity from its very sparse data sampling.

We base our analysis of the HD 27894 system on various models. First we apply a single Keplerian fit just like Moutou et al. (2005) and obtain fit parameters in agreement with those found by these authors (except for a 2 m s^-1 offset in the systemic RV γ; see Table 3). As discussed in Sect. 1, this model has six free parameters (Eqs. (3) and (4)). We find a low-eccentricity fit with a reduced chi-square of ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{8.706}$ (d.o.f. = 14), hence p(χ²) = 2.69 × 10^-19 implying that the fitted model alone cannot at all describe the data if the measurement errors are reliable. At this point, of course, a multi-planet model is not necessarily called for yet, since stellar RV jitter (due to pulsations, convection, and activity) could also be responsible for the excess variability; see the pertinent discussion in Moutou et al. (2005) on the fact that the adopted average RV error of 1.8 m s^-1 does not take stellar variability into account. The single Keplerian fit has an rms of 4.21 m s^-1 and yields a minimum companion mass of msini = 0.645 M_Jup and an orbital period of P = 17.99 d (Table 3).

Next we consider two-planet fits. We start off with the exact 2:1-RCO case which yields values for the orbital period and planetary mass for the outer planet very close to those found for the single-Keplerian model (Table 3) and only a small improvement of chi-square, ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{8.581}$, p(χ²) = 5.92 × 10^-19 (d.o.f. = 14) and rms (4.17 m s^-1). In this model, the inner planet has msini = 8.4 M_Earth. The similarity of the 2:1-RCO fit with the single Keplerian demonstrates the ambiguity of the model selection in this case of sparsely sampled RV data.

The situation changes somewhat when we leave the two period values free but retain circular orbits (eccentricities e₁ = e₂ = 0, hence ω₁ and ω₂ are undefined). Our new model now has seven free parameters, one more than Eq. (1) or (2) since here we are not forcing the 2:1 ratio for the periods. The resulting fit corresponds to a two-planet system with P₁ = 17.9919 d, m₁sini = 0.648 M_Jup, P₂ = 8.276 d, and m₂sini = 0.044 M_Jup = 14 M_Earth (see Table 3 and Fig. 6). The two-planet fit has a substantially smaller ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{5.444}$ (d.o.f. = 13), but still a very small p(χ²) of 2.37 × 10^-19, probably since stellar jitter is dominating the noise behaviour but is not considered in the model. The rms of the fit is 3.20 m s^-1. Applying an F-test to the χ² ratio of the fits, we find a relatively moderate confidence of 83% that this new model is an improvement over either the single Keplerian or the 2:1-RCO models (same confidence value in both cases). If the F-test is applied to the rms’s of the fits, a similar value of 82% is obtained in both cases.

We then study the orbital stability of the possible two-planet system following an approach similar to the one outlined in Trifonov et al. (2014). We perform N-body simulations employing the symplectic integrator SyMBA (Duncan et al. 1998) which is particulary suited to treat close encounters between the planets. The fitting code is described in detail in Tan et al. (2013).

For these simulations we assumed co-planar orbits seen edge-on for both planets, i.e. the orbital inclinations are i₁ = i₂ = 90°. The initial orbital elements are obtained in the Jacobi frame matched to the representation of hierarchical multi-planet systems (see Lee & Peale 2003). Our choice of co-planar edge-on orbits can be motivated as follows. First, if nothing is known about the orbital inclination and random orientation is assumed, 90° is the most probable value. Second, for this inclination value the planet mass is equal to its minimum mass determined from the RV data. While higher mass planets will have a stronger gravitational interaction and can potentially disturb each other more, edge-on orbits are the relevant ones for the question whether there are stable orbits at all. Third, while co-planarity maximizes the interaction between the planets, co-planar orbits are nevertheless thought to be the predominant configuration for planetary systems originating from a circumstellar disk. This leads us to assume them here.

Fig. 7 Orbital evolution of the best-fit system of two planets that started off on circular orbits but with unconstrained periods. Left: 1-year excerpt from the simulation. Right: excerpt for the following 999 yr. Both in the left and right part of the figure, the individual panels show from top to bottom the short-term variation of the semi-major axes a1 and a2, of the eccentricities e1 and e2, and of the difference of the longitudes of periastron ω2 − ω1. The last quantity is not plotted when one of the eccenticities are close to zero (e1 or e2< 0.0001) since then the longitude of periastron is undefined (hence the gaps in the lower left panel). The variation pattern of these parameters shows various periodicites of which the following can be most easily seen here: 15 d and 104 d (all five parameters, left part of figure), and 65 yr (e1 and ω2 − ω1; right part of figure).

Fig. 8 Map of modified chi-square, ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{·}{\mathit{\chi }}^{\mathrm{2}}\mathit{/}{\mathit{\chi }}_{\min }^{\mathrm{2}}$, as a function of the two orbital periods of the two-circular-planets model for HD 27894. All models on this map are stable. Contours are combined confidence levels with respect to the best fit, i.e. the 1, 2, 3 and 4σ regimes in which both period values are found together with a probability of 68.27%, 95.45%, 99.73%, and 99.99%, respectively (while the remainder of the orbital parameters is optimized, i.e. treated as uninteresting). The best fit with ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{13}$ is denoted by an asterisk. It is within the ± 9% region from the exact 2:1 period ratio. The left and right slanted dashed lines delimit this region whereas the middle one marks the exact ratio. Exact 2:1 solutions are found within the 3σ contour (middle). The 1, 2, 3 and 4σ contours correspond to ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{15.30}\mathit{,}\mathrm{19.18}\mathit{,}\mathrm{24.83}$, and 32.33, respectively.

We find the best fit with initially two circular orbits and free period values to be stable over the complete simulation time span of 10 000 yr corresponding to over 2.3 × 10⁵ and 4.6 × 10⁶ orbits of the outer and inner planet, respectively. A 1 yr and a (contiguous) 999 yr excerpt from the simulation are shown in Fig. 7 and reveal very regular oscillations (on several time scales) of the semi-major axes, the eccentricities, and the difference of the longitudes of periastron of the two planets. The amplitudes of these oscillations are small. The semi-major axes of the outer and inner planet vary by about 10^-4 and 2 × 10^-4, repectively. Only very small eccentricity values are reached (0 ≤ e₁ ≤ 0.00045 and 0 ≤ e₂ ≤ 0.013). With a period ratio of P₁/P₂ = 2.17, this system is not exactly in 2:1 MMR but still in the range of period ratios within ± 9% of the MMR found to be frequent by Lissauer et al. (2011) in Kepler data (see Sect. 1 above). It is particularly noteworthy that this period ratio is very close to the value of 2.2 where a significant excess of planet pairs was found by Steffen & Hwang (2015) also in Kepler data (Sect. 1).

In order to check whether a system starting off at period ratios near 2:1 (including the exact ratio) can also be stable, we examine a high density (50 × 50 fits) χ² grid of two-planet circular fits with different initial periods. We have tested each fit for long-term stability. The result is shown in Fig. 8 which shows a modified chi-square defined to be equal to d.o.f. at its minimum by ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{·}{\mathit{\chi }}^{\mathrm{2}}\mathit{/}{\mathit{\chi }}_{\min }^{\mathrm{2}}$. Note that all models on this map are stable, and so is the one starting at the exact 2:1 period ratio which is also consistent with the optimum fit (at ${\mathit{\chi }}_{\min }^{\mathrm{2}}$) at better than the 3σ level; more precisely, it lies within the combined 3σ error regime for the interesting parameters P₁ and P₂ around the optimum fit (cf. Lampton et al. 1976; also Avni 1976).

Finally, we fit the eccentric two-planet model to the data. First we apply the double Keplerian model, but we are unable to constrain the eccentricity of the inner planet (perhaps due to the sparse data sample). The best-fitting double Keplerian has ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{2.522}$ and p(χ²) of 6.91 × 10^-3 (d.o.f. = 9; see Table 3 for the fit parameters). When comparing it with the 6-parameter single Keplerian and 2:1-RCO models, an F-test applied to the χ² ratios shows that in both cases there is a confidence of 99% that the five additional parameters of the double Keplerian provide a significant improvement of the fit. When the F-test is applied to the ratio of the squared rms’s, this confidence is 95% in both cases. However, the best-fitting double Keplerian yields a very large value for the eccentricity of the inner planet of e₂ = 0.51. Not only does this value have an enormeous uncertainty of ± 0.47, it also leads to a highly unstable planetary configuration. In order to exclude such unstable models by imposing stability constraints, we apply a self-consistent two-planet N-body model (Laughlin & Chambers 2001; Tan et al. 2013). We find that, if we constrain the eccentricity of the inner planet to e₂< 0.3, we arrive at a stable orbit again.

Comparing the models we have studied in this subsection we certainly find some preference for a two-planet model for the HD 27894 system even though this is not fully conclusive yet. This is mainly due to the sparseness of the available data. Since the formal best fit with two Keplerians drives the eccentricity of the potential inner planet to values so high that the system becomes unstable, models with moderate eccentricities should be preferred. This includes the model with two circular orbits. HD 27894 is clearly an example for a system where more observations hold the potential to uncover a second planet.

5. Discussion

5.1. Observational strategy for the follow-up

In this paper we have illustrated how RV data with uncertainties can lead observers to apply the wrong model for the description of the data of planet hosting stars. In particular, we have concentrated on the topic of the erroneous application of a single Keplerian one-planet model instead of a 2:1 resonant two-planet model. Sparse sampling as well as insufficient measurement precision can make it difficult to distinguish between the two models. Sparse sampling may just not cover enough of those orbital phases where the two models differ most, and measurement errors may be too large for attempts to resolve the difference between the models.

Naturally, improvement will come from more observations of the candidate stars, and even more so if also the measurement precision can be improved. Observational methods other than RV measurements, such as astrometry or photometric transits can contribute to the selection of the best model for the studied planetary system. McArthur et al. (2014) provide an example for employing constraints from astrometry and photometry to distinguish between or exclude certain models. Mostly though, the follow-up will consist of securing more RV measurements. These can be optimized based on the following recipe:

• Fill phase gaps in the data.

• Preferentially observe those phases of the RV curve where the differences between the models are strongest.

On average, doubling the number of measurements increases the ability to resolve the difference between the models by a factor of $\sqrt{\mathrm{2}}$. With higher instrumental precision this ability is even greater provided that the main limitation is not due to intrinsic effects of the star such as activity, convective motions, and pulsations.

To illustrate the best choice of phases by way of example, take the second panel from the left in the upper row of panels in Fig. 1. If one compares the thick solid line, i.e. the 2:1-RCO two-planet model, with the thin solid line, i.e. the Keplerian one-planet model, the largest differences occur at phases near 0.18, 0.40, 0.60, and 0.82 whereas intermediate phases have smaller discrepancies. Given the usually limited observing time, one should concentrate on those phases of the largest differences.

5.2. Biased eccentricity distribution and consequences for planetary system formantion

Determining the fraction of two-planet systems near the 2:1 resonance among the systems previously classified as eccentric single-planet ones will allow us to find the degree to which the observed eccentricity distribution is biased. It may turn out that solar-system-like architectures with low-eccentricity planets are more frequent than previously thought. If a significant bias can be established, this will certainly have implications for theories of the formation of planetary systems. We note in this context also that, according to several authors, resonances observed in planetary systems provide clues to planetary migration processes and the architecture of those systems (e.g. Pierens et al. 2014; Delisle & Laskar 2014; Beaugé et al. 2006).

5.3. Outer planets as perturbers

One related topic is worth mentioning here, even though we have not studied it in this paper. Rodigas & Hinz (2009) investigated the effect that undetected second planets in wide (outer) orbits have on the determined eccentricity value of supposed single planets in sparse RV data. These authors find that a massive outer planet can lead to published eccentricity values for the detected inner planet that are too large. In this way they have a similar effect as the 2:1 resonant inner planet on which we have concentrated in this study.

6. Conclusions

• 1.

  We have studied the case that in RV measurements a two-planet system with 2:1 resonant circular orbits (2:1-RCO system) is erroneously interpreted as a single-planet system and fitted with a single eccentric Keplerian.

• 2.

  By way of simulation we have determined how, on average, the discrepancy between the true and the false models depends on the eccentricity of the fitted Keplerian.

• 3.

  From the excess scatter of the data around the fitted model, we find that in 74% of 254 stars for which a single eccentric planet was announced a 2:1-RCO system is a possible alternative model that cannot be excluded.

• 4.

  Based on the probability of χ², we find that for 54% of 187 stars the Keplerian fit can be rejected as the sufficient model at the 95% confidence level and for 39% of 187 stars it can be rejected at the 99.9% confidence level.

• 5.

  Therefore, a large number of single-planet systems announced in the literature should be scrutinized with follow-up observations in order to provide the data quality needed to distinguish between the two models in question.

• 6.

  Conceivably, it may then be further established that a substantial fraction of the supposed single eccentric orbits are in fact near-circular systems of two planets near the 2:1 resonance. This finding would imply a bias in the determined eccentricity distribution making the quasi-circular orbits prevailing in the solar system more frequent than presently thought. This would have direct implications on planet formation theory.

• 7.

  Employing the example of HD 27894, we illustrate the model alternatives finding that this star may have a Uranus or Neptune-mass planet in an inner orbit near the 2:1 resonance with its known Jovian planet. Intriguingly, our two-planet models for this star yield a period ratio close to 2.2, a value where an excess of planets pairs is found in data from the Kepler mission.

• 8.

  For this system, we also perform stability simulations showing that there are stable orbits at or near the 2:1 MMR when we assume initially circular orbits for both planets.

• 9.

  If fully Keplerian orbits are allowed for both the known and the hypothesized planet, the HD 27894 system is only stable when the eccentricity of the inner planet is constrained to < 0.3.

• 10.

  We suggest that follow-up observations of candidate systems should naturally concentrate on those phases of the orbit of the known planet with the strongest discrepancies between the two-planet and the single-planet (Keplerian) models.

Acknowledgments

We thank the anonymous referee for useful suggestions. N.M.K. ackowledges support by the European Research Council Advanced Grant HotMol (ERC-2011-AdG 291659). F.R. acknowledges financial support from the Alexander von Humboldt foundation.

References

All Tables

All Figures

Fig. 1 Examples for single-component fits to the RV data of a pair of 2:1-RCO planets. Thick solid lines represent the input data which are shown for amplitude ratios and 2/1 (from left to right) and values of the phase shift Δφ = 0° (upper row of panels) and 90° (lower row of panels). Thin solid lines depict the best-fitting Keplerian orbits with P = 1 (not discernible in the leftmost panels because of their very good match to the thick solid line). Dotted lines indicate best-fitting circular orbits (sine waves) with P = 1. Dashed lines stand for best-fitting circular orbits with P = 0.5. Note that the best-fitting Keplerian orbits with P = 0.5 are identical to the circular orbits with the same period, since the eccentricity vanishes for all fits. Plots for phase shifts of 180° and 270° can be obtained by flipping the 0° and the 90° plots, respectively, around the horizontal axis and then shifting them by 0.5 in phase.

In the text

Fig. 2 Results from our simulations of fitting 2:1-RCO data with various models. Left panel: the rms deviation σmodel between the various types of fit and the 2:1-RCO input data as a function of the amplitude ratio of the input data. Shown is σmodel for the single-component Keplerian fits (“KF”) and for the sine fits (“SF”) pertinent to the circular models. Note that σmodel is expressed in units of the semi-amplitude K of the Keplerian fit with Period P = 1. Middle panel: variation of the eccentricity e of the Keplerian fits with . Right panel: the semi-amplitude K of the Keplerian fits as a function of and normalized to one for . In all panels fits with period P = 1 are plotted with larger symbols than fits with P = 0.5 (see the labels at the curves). Thin dotted lines in all three panels: as a comparison we have included the result from a simulation for sparse data sampling based on the times at which the 20 RV measurements for the host star of HD 27894 were taken by Moutou et al. (2005). These thin dotted lines correspond to the median values for σmodel, e, and K, respectively, obtained from this simulation (Keplerian fits with P = 1 only). Note that for less than ≈ 0.5 they are located near the line for KF (P = 1) for e and K, whereas for σmodel the similarity holds up to .

In the text

Fig. 3 Rms deviation σmodel of the Keplerian fit with P = 1 from the 2:1 resonant input model (solid line with points and labels). This deviation is expressed in units of the RV semi-amplitude K and plotted as a function of the eccentricity e determined by the fit. The labels indicate the amplitude ratio of the input sinusoids. Thin dashed line: amplitude of the first harmonic in the Fourier expansion of the Kepler equation.

In the text

Fig. 4 Rms scatter σKep of the RV solutions around single-planet Keplerian fits for our sample of 254 stars for which the EOD provides these values (dots). The scatter values are plotted as a function of the formal eccentricity of the Keplerian fit and in units of its RV semi-amplitude. 1σ-errors of the eccentricity are represented by grey horizontal bars. In two instances no errors are available; the pertinent data points are marked by small circles around them. The large circle and cross marks the case of HD 27894 that is studied in Sect. 4.3. The solid curve is identical to the solid line in Fig. 3 representing the deviation σmodel between the single-planet Keplerian model and the model for the 2:1 resonant pair of planets. Scatter values above the curve are larger than the expected deviation arising from applying the erroneous model (single Keplerian instead of 2:1 resonant pair) so that both models are possible, whereas scatter values below the curve are too small to assume an erroneous model and therefore favour the single Keplerian.

In the text

Fig. 5 Cumulative histogram of $\log \mathrm{\left[}\mathit{p}\mathrm{\right(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}$ of the Keplerian fits for those 187 stars in our sample for which both the reduced chi-square ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ of the fit and the number of measurements is provided by the EOD. For a substantial fraction of the stars, the chi-square probability is too small to accept the single Keplerian model, e.g. 39% of the stars have $\log \mathrm{\left[}\mathit{p}\mathrm{\right(}{\mathit{\chi }}_{\mathrm{Kep}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}\mathit{<}\mathrm{-}\mathrm{3}$ (see also Table 2).

In the text

Fig. 6 Fit to the HD 27894 RVs with a model consisting of two planets in circular orbits. Note that in this case the exact 2:1 resonance is not forced, but the two period values are both free parameters of the fit. Left: time series together with the best-fitting two-planet circular model (top, thick solid line) and residuals (bottom). For display purposes the x-axis has been broken between BJD – 2 452 950 and BJD – 2 453 267.5 where there is a large data gap. The thin dashed line in the upper panel represents the best-fitting single Keplerian model for comparison. Right: RV data phase folded with the period of the outer planet after subtraction of the signal of the inner planet (top) and phase folded with the period of the inner planet after subtraction of the signal of the outer planet (bottom).

In the text

Fig. 7 Orbital evolution of the best-fit system of two planets that started off on circular orbits but with unconstrained periods. Left: 1-year excerpt from the simulation. Right: excerpt for the following 999 yr. Both in the left and right part of the figure, the individual panels show from top to bottom the short-term variation of the semi-major axes a1 and a2, of the eccentricities e1 and e2, and of the difference of the longitudes of periastron ω2 − ω1. The last quantity is not plotted when one of the eccenticities are close to zero (e1 or e2< 0.0001) since then the longitude of periastron is undefined (hence the gaps in the lower left panel). The variation pattern of these parameters shows various periodicites of which the following can be most easily seen here: 15 d and 104 d (all five parameters, left part of figure), and 65 yr (e1 and ω2 − ω1; right part of figure).

In the text

Fig. 8 Map of modified chi-square, ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{·}{\mathit{\chi }}^{\mathrm{2}}\mathit{/}{\mathit{\chi }}_{\min }^{\mathrm{2}}$, as a function of the two orbital periods of the two-circular-planets model for HD 27894. All models on this map are stable. Contours are combined confidence levels with respect to the best fit, i.e. the 1, 2, 3 and 4σ regimes in which both period values are found together with a probability of 68.27%, 95.45%, 99.73%, and 99.99%, respectively (while the remainder of the orbital parameters is optimized, i.e. treated as uninteresting). The best fit with ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{d}\mathit{.}\mathrm{o}\mathit{.}\mathrm{f}\mathit{.}\mathrm{=}\mathrm{13}$ is denoted by an asterisk. It is within the ± 9% region from the exact 2:1 period ratio. The left and right slanted dashed lines delimit this region whereas the middle one marks the exact ratio. Exact 2:1 solutions are found within the 3σ contour (middle). The 1, 2, 3 and 4σ contours correspond to ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\mathrm{15.30}\mathit{,}\mathrm{19.18}\mathit{,}\mathrm{24.83}$, and 32.33, respectively.

In the text