aa r X i v : . [ a s t r o - ph ] S e p Celestial Mechanics & Dynamical Astronomy manuscript No. (will be inserted by the editor)
Resonances of low orders in the planetary system of HD37124
Roman V. Baluev
Received: April 8, 2008 / Revised: June 21, 2008 / Revised: September 1, 2008 / Accepted: September 10,2008
Abstract
The full set of published radial velocity data (52 measurements from Keck + 58ones from ELODIE + 17 ones from CORALIE) for the star HD37124 is analysed. Twofamilies of dynamically stable high-eccentricity orbital solutions for the planetary systemare found. In the first one, the outer planets c and d are trapped in the 2/1 mean-motionresonance. The second family of solutions corresponds to the 5/2 mean-motion resonancebetween these planets. In both families, the planets are locked in (or close to) an apsidalcorotation resonance. In the case of the 2/1 MMR, it is an asymmetric apsidal corotation(with the difference between the longitudes of periastra Dw ∼ ◦ ), whereas in the case ofthe 5/2 MMR it is a symmetric antialigned one ( Dw = ◦ ).It remains also possible that the two outer planets are not trapped in an orbital resonance.Then their orbital eccentricities should be relatively small (less than, say, 0 .
15) and the ratioof their orbital periods is unlikely to exceed 2 . − . Keywords planetary systems · resonance · stability · periodic orbits · statistical methods For now, the star HD37124 is believed to host three Jovian planets. The innermost planet‘b’ was discovered by Vogt et al. (2000). This planet moves on a low-eccentric orbit with aperiod of P b ≈
150 days. Soon after this, the second planet ‘c’ was discovered independentlyby Udry et al. (2003) and Butler et al. (2003). At that time, its mass and orbital parameters(e.g. period P c ∼ R.V. BaluevSobolev Astronomical Institute, Saint Petersburg State University,Universitetskij prospekt 28, Petrodvorets, Saint Petersburg 198504, RussiaE-mail: [email protected] There is no clear consensus between researchers about notation of planets in the system. The innermostplanet is always denoted as ‘b’, but the notation for the outer pair of planets may vary. We use the samenotation as used in
The Extrasolar Planets Encyclopaedia by J. Schneider, . Namely,we denote the innermost, the outermost and the intermediate planet by the letters ‘b’,‘c’,‘d’, respectively.
Still, the radial velocity (RV) data are insufficient to obtain reliable estimations of pa-rameters of this system directly. Any attempt to obtain a best-fitting RV solution inevitablyleads to dynamically unstable orbital configuration disintegrating after a very short timedue to high eccentricities of two outer planets. To force the fitting algorithm to find a sta-ble configuration, Vogt et al. (2005) fixed the value of e c at 0 .
2. Go´zdziewski et al. (2006,2008) presented a detailed analysis of the Keck RV data involving constraints of dynamicalstability. The stable orbital configurations from these works span a very wide region, withratio of orbital periods of the outer planets, P c / P d , ranging from ∼ . ∼ .
9. We haveto ascertain that, in fact, the Keck RV data only outline a wide region of acceptable orbitalconfigurations, whereas significant constraints on the orbits of these planets are set mainlyby the stability requirement.The aim of the present paper is to describe the set of most likely orbital solutions forthe system of HD37124, based on the analysis of the complete set of RV data published(incorporating Keck, ELODIE, and CORALIE measurements). The structure of the paperis as follows. In Section 2, the RV datasets used in the paper are described. In Section 3,the statistical methods used in the paper are discussed. In Sections 4, 5, and 6, the resultsof the analysis are presented. In Section 7, the dynamical behaviour of the resulting orbitalconfigurations is considered. In Section 8, the hypothesis of existence of an extra planet inthe system is tested.
The most precise publicly available radial velocity data for HD37124 were published in thepaper by Vogt et al. (2005). These 52 measurements were obtained at the Keck telescope andspan about 8 . . . . . ∼ P b ≈
150 days in HD37124 is much longer. The 58 reconstructed ELODIEdata points cover about 7 . . . . . . . J = N j ( j = , , ) RV measurements v ji ( i = , , . . ., N j ) having the‘stated’ RV uncertainties s meas , ji and made at the epochs t ji . These datasets are plotted inFig. 1, top panel.It was demonstrated in (Baluev, 2008b,c) that high-precision RV measurements in planetsearch surveys often suffer from periodic (annual) systematic errors which may originatefrom various sources. Sometimes, these systematic errors may reach the magnitude ∼
10 m/s,especially for data published several years ago, when data reduction algorithms have notbeen debugged to a perfect state yet. It is necessary to account for these systematic errors inour analysis. For this purpose, a simple harmonic model of these errors A cos ( p ( t − t ) / ) R A D I A L VE L O C I T Y [ m / s ] -50 0-10050-50 0-10050-50 0 50 100 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 R A D I A L VE L O C I T Y [ m / s ] Fig. 1
Top panel: RVmeasurements for the starHD37124. Different RVoffsets were assigned to theELODIE, CORALIE, andKeck datasets, in order toseparate them from eachother. Bottom panel: RVcurve for the best stable fitfrom (Go´zdziewski et al.,2008), its RV residuals,and its RV residuals aftercorrection of the best-fittingsinusoidal annual drift in theELODIE data. is adopted below. Here, the semi-amplitude A and the time shift t are the extra free param-eters to be determined from the time series. In the next section we will see that ELODIEdata always show a significant annual drift of radial velocity with large semi-amplitude A ∼
20 m/s. Unfortunately, the number of CORALIE measurements is too small for reli-able modelling of their possible annual errors. Seemingly, it is better to avoid this modellingfor the CORALIE data. The existence of significant annual errors in Keck measurements isuncertain. We will consider different models of the RV curve below, with and without theannual term in the Keck data.One might ask that since the data from ELODIE and CORALIE are less accurate, sufferfrom significant systematic errors, and increase the total time coverage by only 20 − ∼
30 m/s for theplanet d and ∼ −
20 m/s for the planet c. Such arcs of the RV curves are very important forconstraining the corresponding orbital periods, which result in more strict constraining ofthe whole set of parameters. Of course, the Keck data remain the main source of information and drive the fit, but the data from ELODIE are also important, because they can help to ruleout a large fraction of inacceptable fits. This is illustrated by the bottom panel in Fig. 1. Wecan see that the orbital solution from the work (Go´zdziewski et al., 2008) fits satisfactorilyall available RV data in the range 1997–2005, except for the ELODIE measurements madein 1995–1996, which show systematic deviation reaching 50 m/s. This makes the mentionedorbital solution significantly less credible. Other solutions, not being ruled out by the Keckdata, may produce even larger deviations and thus could be easily ruled out even by veryinaccurate measurements. However, data from different observatories have different statis-tical properties and it is necessary to merge them extremely carefully, in order to set correctstatistical weights to different datasets. This problem will be considered in the next section. j th observatoryat time t as m j ( t , p ) = m obs , j ( t , p obs , j ) + m ⋆ ( t , p ⋆ ) , j = , , . . . J . (1)Here, the full vector p of free parameters to be estimated consists of elements of vectors p obs , j ( j = , , . . ., J ) and p ⋆ . The function m obs , j in (1) represent an observatory-specificpart of the measured radial velocity: m obs , j ( t , p obs , j ) = c , j + s j (cid:229) n = A jn cos ( p ( t − t jn ) / P jn ) (2)The constant velocity term c , j and parameters A jn , t jn of possible systematic errors form thevectors p obs , j of unknowns. The quantities P jn — the periods of the systematic errors — areassumed a priori known. In this paper, we will consider only the cases s j = s j = P j being the annual period. The function m ⋆ in (1) is the commonradial velocity term incorporating RV signals due to unseen companions orbiting the star: m ⋆ ( t , p ⋆ ) = r (cid:229) n = c n t n + N (cid:229) n = K n ( cos ( w n + u n ) + e n cos w n ) . (3)The coefficients c n describe possible long-term polynomial (of degree r in general, we willconsider only the cases r =
0, no trend, and r =
1, linear trend, below) trend in the RV data.This trend may be induced by possible distant unseen companions in the system with periodslonger than the time span of the observations. Other terms in (3) represent Keplerian veloc-ities induced by N planets ( N = c n , RV semi-amplitudes K n and the orbital elements l n (the mean longitude at certain fixed epoch), e n (the eccen-tricity), w n (the argument of the periastron), P n (the orbital period) form the vector p ⋆ ofplanetary parameters to be estimated. The quantity u n in (3) is the true anomaly of the n th planet (evidently, it depends on the time and on the parameters l n , P n , e n ).The model (3) does not take into account interactions between the planets. This is ad-missible for the case of HD37124, because after ∼
10 yr of RV observations, the outermostplanets ‘c’ and ‘d’ (which, as we will see below, represent the main source of dynamical ac-tivity in the system) have completed only two and four revolutions around the star. On sucha time scale, the planetary perturbations could be only significant in the case of sufficiently -40-30-20-10 0 10 20 30 40 50 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 N e w t on i an R V - K ep l e r i an R V [ c m / s ] Fig. 2
The difference between the Newto-nian N-body and multi-Keplerian RV modelsfor the fit from (Go´zdziewski et al., 2008).The clear growth of this difference up to ∼
50 cm/s reflects the fact that the orbital el-ements were given for the epoch of the firstKeck observation. If the reference epoch wascloser to the middle of the observation win-dow, the RV difference would be bounded by20 −
30 cm/s over the whole time segment. close approaches between planets on high-eccentricity orbits. However, such configurationsare unrealistic, because close approaches usually represent the source of instability and leadthe system to disintegration. For example, the difference between the best-fitting Keplerianand N-body RV models (see Fig. 2) do not exceed 50 cm/s for the best stable orbital con-figuration from (Go´zdziewski et al., 2008). For stable orbital solutions that are consideredbelow, these deviations do not exceed ∼
30 cm/s and produce only ∼ − relative error inthe RV residuals r.m.s. Bearing this in mind, only the Keplerian RV models are used in thepaper. They allow much faster computations than the N-body ones.The minimum mass of the planet m sin i (here i is the orbital inclination to the sky plane)and its semi-major axis a can be derived as m sin i ≃ ˜ K (cid:18) M ⋆ P p G (cid:19) / = M ˜ KP / M / ⋆ , a ≃ (cid:18) GM ⋆ P p (cid:19) / = A P / M / ⋆ , (4)where M ≈ . · − [ M Jup · M − / ⊙ · m − · day − / · s] and A ≈ . · − [AU · M − / ⊙ · day − / ] are constant factors, G is the gravitational constant, M ⋆ is the mass of the star, and˜ K = K √ − e is the modified semi-amplitude. As it is well known, the inclination i cannotbe constrained using the Keplerian RV model.The errors of the approximate equalities in (4) are about m / M ⋆ ∼ − and are much lessthan the statistical uncertainties of estimations given below. For example, the shift in the RVsemi-amplitude of about m / M ⋆ ∼ − would lead to an error in the RV of only 1 − O ( m / M ⋆ ) . For these cases, we need to provide more decimal digits (excessive with respectto the statistical uncertainties) for the values of the semi-major axes, in order to allow thereader to reproduce the results of long-term numerical integrations discussed below. Forthis purpose we should state the coordinate system in which we refer our estimations. Itwas noted by Lissauer & Rivera (2001) and Lee & Peale (2003) that it is better to interpretthe orbital parameters of the Keplerian model as osculating ones referenced in the Jacobicoordinate system. In the Jacobi coordinates, the osculating orbital period and semi-majoraxis of an k th planet are connected by the relation a k = A P / k M ⋆ + k (cid:229) j = m j ! / . (5) Here we need to substitute the values of the planetary masses themselves, and thus to assumesome values for orbital inclinations. As it will be discussed in Section 5, we assume that theplanetary system is seen edge-on, that is i k = ◦ .3.2 Estimations of the parametersTo analyse the RV data described in Section 2, we need to estimate the so-called RV jitter s ⋆ which increases the full RV uncertainties as s = s ⋆ + s and softens the differ-ences between the weights of observations (cid:181) / s . In the astrophysical part, this RV jitteris inspired by various activity in the star (e.g. Wright, 2005), but often incorporates instru-mental effects as well. As is shown by Baluev (2008c), the effective RV jitter may be quitedifferent for different instruments, even for one and the same star. Therefore, we shouldperform the merging of RV datasets from different observatories very carefully, with correctassignment of statistical weights to these data. To do it, we use here the maximum-likelihoodapproach described in the paper (Baluev, 2008c). This algorithm includes a built-in estima-tion of the effective RV jitter (simultaneous with the estimation of usual parameters), whichallows us not to rely on a low-precision astrophysical estimations of s ⋆ . Moreover, this al-gorithm allows to perform a separate estimation of the effective RV jitters for the datasetsfrom different observatories. This algorithm uses the maximization of the modified log-likelihood function of the N RV observations v ji (their errors are assumed to be uncorrelatedand Gaussian-distributed) which is defined asln ˜ L = − J (cid:229) j = N j (cid:229) i = " ( v ji − m j ( t ji , p )) gs , ji + ln s full , ji − N ln √ p s ⋆ , p −→ max . (6)Here, s , ji = s ⋆, j + s , ji ( j = , , . . ., J ; i = , , . . . N j ) and the correction divisor g = − d / N with d being the number of degrees of freedom in our RV model (i.e., the numberof free parameters, d = dim p ). The divisor g allows to perform a ‘preventive’ reduction ofthe statistical bias in the estimations of the RV jitter. Since we aim to use the new objectivefunction ˜ L , instead of the usually used c one, we need to introduce a new measure of thegoodness-of-fit, which would be based on ˜ L . In accordance with Baluev (2008c), to assessthe quality of a given orbital fit, the following goodness-of-fit measure is used below:˜ l = ˜ L − / N e − . / √ p ≈ . L − / N . (7)This function is measured in the same units as radial velocity (i.e., in m/s). It characterisesthe overall scatter of RV measurements around the model. However, to allow a comparisonwith previous works on HD37124, the traditional r.m.s. goodness-of-fit measure is usedbelow as well.3.3 Assessing the reliability of orbital fitsIt is not enough to find an RV fit with a small scatter of residuals. To interpret the resultingestimations, we need to assess their reliability. It is shown by Beaug´e et al. (2008) that orbitalfits of multi-planetary systems in a mean-motion resonance (hereafter MMR) may be highlyunreliable, though formal uncertainties of estimations may be apparently small. In thesecases, the shape of the likelihood function may be complicated and may possess multiplecomparable local maxima. Often this shape is model-dependent: addition of extra model components leads to qualitative changes of the set of likelihood maxima. Often, every suchlocal maximum provides a good fit of the RV curve but nevertheless is far from the actualorbital configuration. Such situation indicates one of the following items:1. The adopted model is imperfect. Some extra terms were not included or the terms in-cluded are wrong.2. The data are imperfect. The errors may have a non-Gaussian distribution, they may becorrelated or they may incorporate some extra time-variable systematic part. Also, thedata may cover too small time base or simply the number of observations is too small.Formally, we could use one of the local maxima (e.g. the global maximum) to construct anestimation of the parameters of the system. However, such estimation would appear stronglybiased (usually to higher eccentricities) and its formal uncertainties would strongly under-estimate real errors of parameters. Therefore, in addition to the formal goodness-of-fit mea-sure, we need some indicator of statistical reliability of our fits.Beaug´e et al. (2008) performed (for the system of HD82943) several fits with truncatedRV datasets to explore the sensitivity of current orbital fits to future RV measurements.Unfortunately, this approach require too time-consuming computations. Here we need somesimple and rapid (though perhaps quite rough) test of the ‘statistical health’ of our orbitalfits. For this goal, we use below the following approach. Recall (Lehman, 1983, § N → ¥ ) approximation to the variance-covariance matrix of the estimationsof p is calculated as the inverse of the Fisher information matrix Q having elements Q ab = J (cid:229) j = N j (cid:229) i = s , ji ¶m j ¶ p a (cid:12)(cid:12)(cid:12)(cid:12) t = t ji ¶m j ¶ p b (cid:12)(cid:12)(cid:12)(cid:12) t = t ji . (8)When we deal with a well-conditioned situation, the likelihood function can be quadraticallyapproximated in the vicinity of its maximum using the quadratic term (cid:181) d p T Q d p only. Thisis the asymptotic behaviour of the likelihood function for large time series (i.e., N → ¥ ).The mentioned quadratic term approximates the multidimensional graph of the likelihoodfunction by a paraboloidal hypersurface. Other terms in the expansion of the likelihood func-tion are insignificant and do not distort this shape essentially. However, when the number ofobservations is not sufficient for reliable determination of the parameters of the model, thematrix Q become ill-conditioned, and the extra terms easily produce distortions leading tomultiple local maxima of the likelihood function.To assess the reliability of a particular fit we could calculate the condition number ofthe information matrix Q (i.e. the ratio of the biggest eigenvalue and the smallest one).The larger is this number, the lower is the reliability of the corresponding orbital fit. Thiscondition number should be compared with the number of observations, because high-orderterms in the expansion of the likelihood function tend to zero when N grows. To demonstratethis, let us write down the full expression of the Hessian matrix of the usual c function(speaking more accurately, of the weighted average of squared residuals): H ab = − Q ab + J (cid:229) j = N j (cid:229) i = m j − v ji s , ji ¶ m j ¶ p a ¶ p b (cid:12)(cid:12)(cid:12)(cid:12) t = t ji . (9)In this expression, the first term determines the asymptotic behaviour of the likelihood func-tion in the vicinity of its maximum. The second term (summation) reflects the non-linearityof the RV model and contains the residuals ( v − m ) , which average decreases as O ( / √ N ) when N grows. This implies that the perturbing term in (9) grows according to the law O ( √ N ) , whereas the Fisher matrix grows according to O ( N ) . Therefore, the relative mag-nitude of the second term in (9) is about 1 / √ N . Let us now transform Q to its diagonal formwith eigenvalues E . . . E d on the diagonal (we assume that E is the biggest eigenvalue and E d is the smallest one). We may expect that the typical magnitude of the elements of theFisher matrix is of the order of √ E E d and, therefore, the typical magnitude of the secondterm in (9) is about √ E E d / √ N . The full matrix H is well approximated by its asymptoticrepresentation − Q , if relative deviations of the corresponding eigenvalues are small. If thisis so, the topology of the graph of the likelihood function should not be very sensitive tothe perturbing terms. We can expect that absolute deviations of the eigenvalues should havesimilar magnitude about p E E d / N . But the relative deviation of the smallest eigenvalue E d is as large as p E / E d / √ N . Therefore, the condition number E / E d of the Fisher informa-tion matrix should not exceed the number of observations N for a statistically robust orbitalfit. This is simply a specification of the general fact, known from the numerical analysis,that the condition number is a measure of sensitivity of a matrix to small perturbations.There is a small clause. The elements of Q are measured in different physical units andit would be incorrect to use the condition number of Q itself (since it would depend on thechoice of the measurement units). Instead, we may use the condition number C of the scaledinformation matrix ˜ Q having elements ˜ Q i j = Q i j / p Q ii Q j j .Unfortunately, the indicator C depends on the epoch for which we obtain the fit, espe-cially when the RV model incorporates long-term trends. To obtain an informative estima-tion of C , we should choose the reference epoch near the middle of the time series span(or, when we merge datasets with different characteristics, near the weighted average ofthe timings t i j ). In this position, the value of C will be close to its minimum and, simul-taneously, the set of parameters will possess the best statistical properties (e.g., the mutualcorrelations of different parameters will be minimized). Even if absolute value of the quan-tity C represents a too rough indicator of the fit robustness, this indicator seems to workquite well for the goal of intercomparison of different orbital fits for the same planetarysystem. To illustrate the use of the indicator C , we give here its (minimum) values for sev-eral planetary systems with well-determined orbits: C =
10 with N =
409 (51Peg), C = . N =
109 (70Vir), C = . N =
203 (14Her). Less perfect cases: C =
17 with N =
75 (HD69830), C =
120 with N =
487 (55Cnc). For the system of HD82943, we ob-tain C =
200 with N =
165 (three-planet model with the eccentricity of the outermost planetfixed at zero).
From now on, we consider three main models of the RV data for HD37124. All modelsincorporate three common Keplerian terms as in the eq. (3), the constant velocity terms(separate for the three datasets), and the annual term for the ELODIE dataset. The first RVmodel (I) does not incorporate anything else and thus has d =
20 degrees of freedom. Thesecond one (II) incorporates also an annual term for the Keck dataset and has d =
22. Thethird one (III) incorporates the same terms as (I) plus a linear trend (common for all the threedatasets) and thus has d =
21. For all these models, the ratio d / N ≈ .
16 means that we areleft with only ≈ d ) free parameters x , y (say the O R B I T A L PE R I O D P _d [ da ys ] O R B I T A L PE R I O D P _d [ da ys ] Fig. 3
Contour maps of the likelihood goodness-of-fit statistic min ˜ l for RV fits of HD37124. The maps areplotted in the plane of orbital periods P c and P d (measured in days). For each point in these panels, thelikelihood goodness-of-fit function ˜ l was minimized over the rest of free parameters. Thick straight line markthe position of the 2/1 MMR. The level values of the plotted function are shown below the respective panels.The panels in different columns correspond to the different RV models (model I, model II, model III, fromleft to right). The top raw of panels corresponds to the fits based on the whole accessible data array, whereasthe bottom raw corresponds to the fits based on the Keck data only. Note the differences in period rangesshown in top and bottom graphs. orbital periods P c and P d of the outer planets). Then we consider the function ˜ l ′ ( x , y ) = min ˜ l , where the minimization of the goodness-of-fit function ˜ l is performed over the restof free parameters. This means that for any manually assigned values of x , y we obtaincorresponding best-fitting values of other parameters and find corresponding goodness-of-fit measure ˜ l . Further, the resulting function of two variables can be visualised on a two-dimensional grid.Fig. 3 shows such plot in the plane of orbital periods P c and P d . We can see that thelikelihood function constructed from the full available dataset (top panels of Fig. 3) has twomain maxima with comparable values of ˜ l . The first one is centred on P c ≈ P d ≈
870 days, and the second one on P c ≈ P d ≈
900 days. We can see that thelatter solution is close to the 2/1 MMR of the outer planets. Such orbital configurations areremarkable because only low-order MMRs can prevent planets on high-eccentricity orbitsfrom close approaches and hence can make the whole system stable. This ‘double-headed’shape of the likelihood function looks rather stable with respect to the choice of the RVmodel. From now on, we consider mainly the two mentioned families of solutions: the firstone corresponds to the resonance 2/1 between planets ‘c’ and ‘d’, the second one is outside Table 1
Best-fitting orbital solutions for the planetary system around HD 37124.parameter IA, 2/1 IB IIA, 2/1 IIB IIIA, 2/1 IIIBplanet b P [days] 154 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) ˜ K [m/s] 28 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) l [ o ] 117 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 131 ( ) ( ) ( ) ( ) ( ) ( ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) planet d P [days] 912 . ( . ) . ( . ) . ( . ) . ( . ) ( ) . ( . ) ˜ K [m/s] 13 . ( . ) . ( . ) . ( . ) . ( ) . ( . ) . ( . ) l [ o ] 290 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 283 ( ) ( ) ( ) ( ) ( ) ( ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 1 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) planet c P [days] 1782 ( ) ( ) ( ) ( ) ( ) ( ) ˜ K [m/s] 14 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) l [ o ] 319 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 256 . ( . ) . ( . ) . ( . ) . ( . ) ( ) . ( . ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 2 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) c [ ms · yr ] - - - - 0 . ( ) . ( ) ELODIE dataset c [m/s] 65 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) A [m/s] 17 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) t [days] 167 ( ) ( ) ( ) − ( ) ( ) ( ) s ⋆ [m/s] 6 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 11 .
85 10 .
96 12 .
23 11 .
21 11 .
79 10 . c [m/s] 4 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) s ⋆ [m/s] 11 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 17 .
45 16 .
93 17 .
62 16 .
73 17 .
68 17 . c [m/s] 7 . ( ) . ( ) . ( . ) . ( . ) . ( ) . ( ) A [m/s] - - 4 . ( . ) . ( . ) - - t [days] - - − ( ) ( ) - - s ⋆ [m/s] 2 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) r.m.s. [m/s] 3 .
704 3 .
713 3 .
044 3 .
381 3 .
566 3 . d
20 22 21˜ l [m/s] 8 .
269 8 .
025 7 .
816 7 .
888 8 .
180 7 . C
213 54 190 74 145 81The values of c for ELODIE and CORALIE are given relatively to their first measuments. The uncertaintiesof the estimations are given in the brackets (e.g., 0 . ( ) means 0 . ± .
10, and 30 . ( . ) means 30 . ± . l and time shift parameters t are given for the epoch JD m sin i and of the semi-major axes a incorporate the 10% uncertaintyof the stellar mass. The estimations of the effective RV jitters s ⋆ incorporate an analytic correction of thestatistical bias as discussed in (Baluev, 2008c).1 of this resonance (but may cover some other MMRs of relatively low, e.g. 7/3 and 5/2).Hereafter, we use the notation ‘A’ for the first family and ‘B’ for the second one.The full sets of estimated parameters for the three RV models are shown in Table 1.The respective minimum values of ˜ l depend on the model adopted. For the models I and III,the ‘B’ solution provides formally better fit to the RV data in comparison with the ‘A’ one,but for the model II the corresponding values of ˜ l are similar. Either model II or model IIIprovide some improvement in the goodness-of-fit ˜ l , with respect to the model I.Let us note that the similar structure of the likelihood function can be also seen in thegraphs constructed in the similar way for the RV model I, but with the use of the Keck dataonly (left-bottom panel in Fig. 3). However, the shape of the likelihood surface is more com-plicated in this case: the broad (much broader than in the top pictures) peak correspondingto the B-family is actually ‘double-headed’ itself (i.e., it consists of two close peaks having P d ≈
850 days and P d ≈
820 days). Moreover, the shape of the likelihood function con-structed with the use of the Keck data only is severely dependent on the adopted RV model.For the model II, the B-family peak is clearly split into two distinct and very distorted peaks(one with P c ≈ P d ≈ P c ≈ P d less than800 days, see middle-bottom panel in Fig. 3). For the model III, all the former peak aremerged into a single very broad peak (formally centred not far from the the 2/1 MMR),which is continued to infinite P c (right-bottom panel in Fig. 3). Such behaviour, the furthersplitting of the local likelihood maxima and their severe sensitivity to the RV model, indi-cates that the analysis of the Keck data alone would yield significantly less reliable resultsthan the joint analysis of all available data.However, even with the use of the full RV dataset, almost all of the best fits possesslarge values of the condition number C , especially for the ‘A’ solution. This means that weshould treat our results with a more care. In fact, no one of the best (in the sense of thegoodness-of-fit measure) fits from Table 1 can be accepted. The values of the eccentricity e c (and often those of the e d as well) are large and lead to a very soon disintegration of thecorresponding orbital configurations. For comparison, Table 2 contains the estimations ofparameters for the system with the eccentricity e c or both the eccentricities e c , e d fixed atzero. These fits have much smaller values of C , though worse goodness-of-fit. Numericalintegration showed their dynamical stability and regular evolution on the time scales of (atleast) 10 yr, except for the fit II ′ A which showed some signs of chaoticity at the time scaleof ∼ yr, evidently due to a large e d = . e c decreases.From the interplay of the indicators C and ˜ l described above, we derive the followingconclusion. Although the orbital fits from Table 1 show relatively small scatter of the resid-uals, this small scattering is in fact fictitious, as indicated by the corresponding values of C .The number of RV data points and their temporal coverage still cannot constrain the full setin the model parameters reliably. The full RV model is ‘overloaded’. Injecting some kindof a priori information (e.g., fixing the eccentricities at low values) allows to overcome theobstacle of the statistical ill-determinacy. The resulting orbital fits possess better statisticalreliability, but by the cost of some increase of the RV residuals scatter. Nevertheless, thisincrease is necessary to obtain a physically realistic orbital configuration.However, we cannot rule out the possibility that the actual orbits of the planets ‘c’,‘d’are far from circular. To find more realistic orbital configurations than those from Table 1,but corresponding to eccentric orbits, we need to account for more subtle requirements ofthe dynamical stability in our analysis. Table 2
Low-eccentricity orbital solutions for the planetary system around HD 37124.parameter I ′ A, 2/1 I ′′ A, 2/1 II ′ A, 2/1 II ′′ A, 2/1 III ′ A, 2/1 III ′′ A, 2/1planet b P [days] 154 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) ˜ K [m/s] 28 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) l [ o ] 117 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 132 ( ) ( ) ( ) ( ) ( ) ( ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) P [days] 908 ( ) . ( ) . . ( . ) ( ) . ( ) . . ( . ) ˜ K [m/s] 15 . ( . ) . ( ) . ( . ) . ( ) . ( ) . ( ) l [ o ] 322 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) ( fixed ) . ( ) ( fixed ) . ( ) ( fixed ) w [ o ] 349 ( ) - 1 . ( . ) - 358 ( ) - m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 1 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) P [days] 1810 ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . K [m/s] 13 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( ) l [ o ] 310 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 2 . ( ) . ( )
61 2 . ( ) . ( )
19 2 . ( ) . ( ) c [ ms · yr ] - - - - 1 . ( ) . ( ) ELODIE dataset c [m/s] 64 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) A [m/s] 18 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) t [days] 167 ( ) ( ) ( ) ( ) ( ) ( ) s ⋆ [m/s] 6 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 11 .
95 12 .
33 12 .
53 12 .
45 12 .
08 12 . c [m/s] 5 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) s ⋆ [m/s] 11 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 17 .
42 16 .
76 17 .
68 16 .
53 17 .
67 17 . c [m/s] 7 . ( ) . ( ) . ( ) . ( . ) . ( ) . ( ) A [m/s] - - 4 . ( . ) . ( . ) - - t [days] - - 134 ( ) ( ) - - s ⋆ [m/s] 3 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) r.m.s. [m/s] 4 .
218 4 .
319 3 .
436 4 .
193 3 .
394 4 . d
18 16 20 18 19 17˜ l [m/s] 8 .
710 8 .
776 8 .
280 8 .
769 8 .
065 8 . C
34 29 40 32 36 31The same notes as in Table 1 to be applied here. In these fits, the eccentricity e c or both the eccentricities e c and e d were fixed at zero. All these configurations correspond to the 2/1 MMR. The estimations of thesemi-major axes and orbital periods are given with one or two excessive decimal digits (shown after two-digit uncertainties enclosed in brackets) to allow an unambigious reproducing of the long-term dynamics ofthese configurations (see Section 3.1). For example, 1810 ( ) . . ±
41 and 2 . ( )
19 means2 . ± . To obtain more realistic stable orbital configurations for this planetary system, we continueto use the method of planar plots of partially minimized goodness-of-fit statistic ˜ l . But nowwe examine orbital solution from a two-dimensional grid more carefully: for each solution,we perform a numerical integration in order to rule out rapidly disintegrating configurations.To perform such integrations, we need to know true masses of planets in the system. As itcan be seen from (4), they depend on the mass of the star M ⋆ and on the orbital inclina-tions. Following Vogt et al. (2005), we adopt M ⋆ = . M ⊙ with an uncertainty of 10%.Unfortunately, using the Keplerian RV model we can estimate only the minimum masses m sin i where the inclination i remains unknown. Until the gravitational interactions betweenplanets in the system are directly observed in the RV curve, the best thing that we can do isto assume a priori that the orbits are coplanar with i = ◦ . If the actual inclination is lessthan 90 ◦ , the true masses of the interacting planets become larger and the stability region ofthe system become more narrow than for the edge-on configurations. The same effect is ex-pected from non-zero mutual inclinations of the orbits, due to the well-known phenomenonof the e – i coupling.Since much troubles in obtaining a realistic orbital configuration of the system are dueto the large eccentricity of the outermost planet, let us firstly consider the plane of eccentricvariables ( e c cos w c , e c sin w c ) . The corresponding maps are plotted in Fig. 4. We can seeclearly the sophisticated shape of the likelihood surface: among the ‘A’ families of solu-tions, no one possess a single maximum. Instead, we can see 2-3 local maxima; all of themcorrespond to large values of e c . The ‘B’ family shows single maximum, but again at large e c . No one of these local maxima corresponds to a stable configuration. Stable solutions oc-cupy only regions of small or moderate e c . It is important to note that the small-eccentricitysolutions correspond to the ‘A’ configuration only: when e c decreases, a B-type solutionapproach the 2/1 resonance and softly turns into an A-type one.As it can be seen from Fig. 4, we can easily find stable solutions from the ‘A’ family.Such solutions can possess values of ˜ l as small as ≈ . ≈ . ≈ . l in the ‘A’ layer, located in the third quadrantof the coordinate plane. On contrary, we face much difficulties with obtaining a stable con-figuration from the family ‘B’. The main reason for almost all ‘B’ solutions to be unstableis that the best-fitting period ratio P c / P d ≈ . − . l is plotted in the plane ( P c , e c ) . We can see that the fits with low e c and with P c fixed far from the 2/1 resonance, possess uncomfortably large values of˜ l . Fig 6 illustrates this more clearly. In this figure, we plot the function ˜ l , minimized sothat the eccentricity e c was fixed at a safe value of 0 .
15 and the period ratio P c / P d wasfixed at the values marked on the abscissas. Since the value of the eccentricity was fixedat a small value, the effects of the RV model non-linearity are significantly suppressed (itfollows from relatively small values of C in Table 2), and thus we may try to find someconfidence intervals for the ratio P c / P d using the standard asymptotic ( N → ¥ ) theory ofpoint estimations. We can see that the values P c / P d > . P c / P d < .
87 are beyondthe 99% confidence interval, when we use the RV model I. For the other RV models, thisconfidence interval becomes even smaller: P c / P d ∈ [ . , . ] for the model II and P c / P d ∈ [ . , . ] for the model III. When the eccentricity e c decreases, these confidence regionsare shrinking around P c / P d = -0.6 -0.4 -0.2 0 0.2 0.4 0.6-0.6-0.4-0.2 0 0.2 0.4 0.6 e_ c s i n ( o m ega_ c ) e_ c s i n ( o m ega_ c ) e_ c s i n ( o m ega_ c ) e_ c s i n ( o m ega_ c ) e_ c s i n ( o m ega_ c ) e_ c s i n ( o m ega_ c ) Fig. 4
Contour maps of the likelihood goodness-of-fit statistic for RV fits of HD37124. These maps areplotted on the plane ( e c cos w c , e c sin w c ) in the same way as in Fig. 3. Solutions which correspond to orbitalconfigurations disintegrating in less than 30000 years are marked by bold grey dots, other solutions aremarked by fine (green in the electronic version of the paper) dots. These dots are arranged according to thepolar coordinate system. Plots in the left column correspond to the ‘A’ solutions close to the 2/1 resonance(with no more than 5% relative deviation of the period ratio), plots in the right column correspond to the ‘B’solutions. The white regions (where e c < .
7) mark the points for which the fitting algorithm could not find asolution from the corresponding family and switched to another one (that appeared significantly more likely).The ‘A’ and ‘B’ families of solutions overlap in the region of large e c without possibility of any smooth seam,but they seem to be sewed smoothly in the region of small e c . The top, middle, and bottom pairs of panelscorrespond to the RV models I, II, and III, respectively.5 E CC E N T R I C I T Y e_ c E CC E N T R I C I T Y e_ c E CC E N T R I C I T Y e_ c Fig. 5
Contour maps ofthe likelihood goodness-of-fit statistic for RV fits ofHD37124. These maps areplotted on the plane ( P c , e c ) in the same way as in Figs. 3and 4. Three panels corre-spond to the models I, II,and III (from top to bot-tom). In each of these pan-els, the A and B families ofsolutions were merged in asingle plot. The bold nearlyvertical lines mark the so-lutions having ratio of thebest fitting periods P c and P d close to the 2/1, 7/3, 5/2,and 8/3 commensurabilities(lines from left to right). L I KE L I H OO D GOO DN ESS - O F - F I T [ m / s ] RATIO P_c / P_d 80% 95% 99% 8.3 8.35 8.4 8.45 8.5 8.55 8.6 8.65 8.7 1.8 1.9 2 2.1 2.2 2.3RATIO P_c / P_d 80% 95% 99% 8.05 8.1 8.15 8.2 8.25 8.3 8.35 8.4 8.45 8.5 1.9 2 2.1 2.2 2.3 2.4RATIO P_c / P_d 80% 95% 99%
Fig. 6
Graphs of the likelihood goodness-of-fit function ˜ l , which was minimized so that the eccentricity e c was fixed at 0 .
15 and the ratio of orbital periods P c / P d was fixed at the values shown on the abscissas. Threepanels correspond to the RV model I, II, and III (from left to right). The bold horizontal lines shows the levelsof min ˜ l yielding the 80%, 95%, and 99% asymptotic confidence intervals for P c / P d . These levels correspondto the values of the likelihood ratio statistic, which provide the asymptotic false alarm probabilities as smallas 20%, 5%, and 1% (see Baluev, 2008c).6 However, we can note a promising region of high-eccentricity solutions near the reso-nance 5 / P c fixed at 2170 days and e c fixed at 0 . years. The evolution of this orbital configuration is a large-amplitude oscillation around anantialigned apsidal corotation, Therefore, this solution belongs to the class of orbital config-urations in the 5/2 MMR found in (Go´zdziewski et al., 2006) using only the Keck data.It seems that, due to the non-linearity of the RV models coupled with lack of the dataand insufficient time coverage, none of the local minima of ˜ l lies near the real orbital config-uration of this system. Probably, these multiple local minima are only the fictitious ‘ripples’produced by the lack of the observations. All these ‘ripples’ are located deeply in the zone ofdynamical instability, and therefore cannot be close to the real configuration of the system.It is possible to find strictly the best-fitting orbital solution(s) simultaneously satisfying thestability requirement. Evidently, such solutions would be attracted by one of the ‘ripples’,and therefore would be close to the boundary of the domain of system stability. Therefore,we would have to deal with large difficulties concerning the very complicated structureof the parameter space near such boundaries (see, e.g., Go´zdziewski et al., 2008), whichprobably represents the Arnold web (see, e.g, Froeschl´e et al., 2006). Near such boundaries,the dynamics of the system is very sensitive to small changes of parameters, the stabilitymap of the parametric space is strongly irregular, and hence the resolution of the parame-ter space should be chosen dense enough. This requires very time-consuming calculationsfor checking the stability of probe orbital configurations in this region. On contrary, it doesnot look likely that real planetary systems can be found in such extremely dynamically ac-tive regions: it would be rather difficult to explain how the system could migrate (withoutdisintegration) to such a state through the dense web of the instability threads and why itstopped in a thin island of stability instead of moving further to dynamically unstable con-figurations. Therefore, the reliability of such ‘hardly-stable’ solutions would be too low tojustify the associated time-consuming calculations. In addition, considering the boundary ofthe stability domain, it is rather difficult to understand the physical mechanism stabilizingthe configurations in the given domain.In this paper, following Hadjidemetriou (2006), we will pay more attention to the cen-tres of the stability domains which point out families of orbital configurations having some‘stock of stability’. For having a clear picture of possible dynamical regimes of a planetarysystem, it is the position of the centre of the stability domain which should be located andfor which we should know possible uncertainties. Considering the kernel of a stability do-main, we avoid dealing with the sophisticated dynamical structure near its boundary, whichis hardly able to carry much information about the dynamics of the real system. When anorbital solution has some stock of stability, the dynamics of the corresponding planetaryconfiguration is much more regular and much less sensitive to small changes of parameters.To obtain such orbital solutions, we will try to decrease the dimension of the problem (i.e.,the number of degrees of freedom d ) using certain a priori information about the stabilityof resonance planetary systems. Bearing in mind the results by Ji et al. (2003); Hadjidemetriou (2006, 2008); Voyatzis & Hadjidemetriou(2006), let us recall that regular stable motions on high-eccentricity orbits with small periodratio are only possible if the planets are trapped in a MMR and simultaneously are close to (or, at least, not far from) an apsidal corotaion resonance. The details of the theory ofapsidal corotation resonances, along with necessary formulae and further references can befound, for instance, in (Beaug´e et al., 2003). For a brief summary, let us consider two planetstrapped in the p / q MMR, i.e. having the ratio of orbital periods P / P ≈ p / q with p > q .We can write down the resonant angles s = p l − q l p − q − w , s = p l − q l p − q − w (10)and the canonically conjugated action variables I = L (cid:18) − q − e (cid:19) , I = L (cid:18) − q − e (cid:19) , (11)where l i are the mean longitudes of the planets and L i ≃ m i √ a i are the Delaunay actionvariables. After averaging the Hamiltonian H of the system over the fast variables (i.e, overthe mean longitudes l i ) keeping the resting slow ones, the resulting averaged Hamiltonian h H i depends on the canonical variables s i , I i , and (as on parameters) on the masses m i of theplanets. Evidently, this averaging accounts properly for the orbital resonance. The averagedequations of motion are then given by dI i dt = − ¶ h H i ¶ s i , ds i dt = ¶ h H i ¶ I i (12)Suppose that some values s ∗ i , I ∗ i determine the position of an extremum of h H i . We can easilysee that every such extremum provides a stationary solution s i ≡ s ∗ i , I i ≡ I ∗ i of the averagedsystem (12). Such stationary solution is often called ‘apsidal corotation resonance’ (here-after ACR). If the initial state of the planetary system slightly deviates from an exact ACR,the motion is a stable oscillation around the exact stationary solution, because the planetsare prevented from close approaches. If the orbits of planets are highly eccentric and are farfrom stationary solutions of the averaged Hamiltonian equations, the motion is, most proba-bly, highly chaotic and unstable: the secular drift of resonant angles (10) leads the planets totoo close approaches destabilizing the system. Stable solutions with one or both the resonantangles circulating are also possible in some cases; however, for high-eccentricity configura-tions, the ACRs mark centres of dynamical stability (see e.g. Hadjidemetriou, 2006, 2008).To obtain nominal orbital configurations of the system, we require from the resonantplanets ‘c’ and ‘d’ to be locked in an exact ACR (while neglecting the influence of the inner-most planet ‘b’). This can be justified not only by the stability considerations. Beaug´e et al.(2006) showed that adiabatic dissipative perturbations (e.g., interaction with a protoplane-tary disk) can cause planet pairs in a MMR to be captured in an ACR as well.The requirement of the ACR lock implies four algebraic equations: ¶ h H i / ¶ s c , d = ¶ h H i / ¶ I c , d =
0. We neglect here the gravitational influence of the innermost planet ‘b’: it isseemingly non-resonant with the outer planets and probably should not affect their resonantdynamics much. The four equations mentioned above put certain constraints on the full set offree parameters to be estimated from RV data and decrease the number of degrees of freedomby four. It is very important, because this decreasing makes the problem significantly betterdetermined: we have about 8 observations per a degree of freedom instead of about 6.More detailed description of the algorithm used to obtain such ACR fits, is given in Ap-pendix A. The resulting fits of ACR solutions are given in Table 3. It is important that sincethe position of the ACR depends on the planetary mass ratio only and is almost independentof the individual planet masses, the solutions from Table 3 are almost independent of the Table 3
ACR fits for the planetary system around HD 37124.parameter I c A, 2/1 I c B, 5/2 II c A, 2/1 II c B, 5/2 III c A, 2/1 III c B, 5/2planet b P [days] 154 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) ˜ K [m/s] 29 . ( ) . ( . ) . ( ) . ( . ) . ( ) . ( . ) l [ o ] 119 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 130 ( ) ( ) ( ) ( ) ( ) ( ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) P [days] 907 . ( . ) . ( . ) . ( . ) . ( . ) ( ) . . ( . ) ˜ K [m/s] 16 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) l [ o ] 315 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 320 . ( . ) . ( . ) . ( . ) . ( . ) ( ) . ( . ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 1 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) P [days] 1815 ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . K [m/s] 12 . ( . ) . ( ) . ( . ) . ( ) . ( . ) . ( ) l [ o ] 311 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) e . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) w [ o ] 267 ( ) . ( . ) ( ) . ( . ) ( ) . ( . ) m sin i [ M Jup ] 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) a [AU] 2 . ( ) . ( )
69 2 . ( ) . ( )
66 2 . ( ) . ( ) c [ ms · yr ] - - - - 0 . ( ) . ( ) ELODIE dataset c [m/s] 64 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) A [m/s] 18 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) t [days] 169 ( ) ( ) ( ) ( ) ( ) ( ) s ⋆ [m/s] 7 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 12 .
20 12 .
60 12 .
85 12 .
51 12 .
15 12 . c [m/s] 4 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) s ⋆ [m/s] 11 . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) . ( . ) r.m.s. [m/s] 17 .
91 15 .
91 17 .
43 15 .
82 18 .
44 16 . c [m/s] 7 . ( ) . ( ) . ( . ) . ( . ) . ( ) . ( ) A [m/s] - - 3 . ( . ) . ( . ) - - t [days] - - 131 ( ) - - - s ⋆ [m/s] 3 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) r.m.s. [m/s] 4 .
322 4 .
108 3 .
762 4 .
148 3 .
910 4 . d
16 18 17˜ l [m/s] 8 .
851 8 .
616 8 .
601 8 .
693 8 .
530 8 . C
39 56 38 64 48 79The same notes as in Table 1 to be applied here. These orbital elements were obtained for the Jacobi coordi-nate system with appropriate masses assigned to the reference barycentres. See Appendix A for more detailsconcerning the algorithm used to obtain these fits. The semi-major axes and some orbital periods are givenwith excessive decimal digits, as in Table 2.9 o -180 o -120 o -60 o o o o o -180 o -120 o -60 o o o o
180 RESONANT ANGLE s_d R ES O N A N T A N G L E s _ c o -180 o -120 o -60 o o o o o -180 o -120 o -60 o o o o
180 RESONANT ANGLE s_d R ES O N A N T A N G L E s _ c Fig. 7
Top panels: orbits of the HD37124 system for the ACR fits (model I only). Left panel is for the 2/1MMR; right panel is for the 5/2 MMR. Straight solid lines mark the conjunction positions. Broken lines markthe lines of apses. Middle panels: contour maps of the averaged Hamiltonian h H cd i in the planes of resonantangles s d , s c , for the same MMRs and eccentricities as in top panels. Crosses mark the actual ACR positionscorresponding to local maxima of h H cd i . assumptions about the orbital inclination of the system: ( m c sin i ) / ( m d sin i ) = m c / m d . Weonly need to adjust the ratio of orbital periods by the quantity O ( m / M ⋆ ) (see formula (22)in Appendix A), but this adjustment have insignificant effect on the RV fit quality (untilwe consider inclinations as small as a few degrees). This invariance with respect to orbitalinclination represents an extra advantage of the use of the ACR fitting procedure.We can see that in the case of the 2/1 resonance, the best fitting ACR is asymmetricwith difference between the longitudes of periastra about 60 ◦ , whereas in the case of the5/2 resonance, the corresponding ACR is symmetric and antialigned (i.e., w c − w d = ◦ ),see Fig. 7. The corresponding RV curves fit all available RV data satisfactorily, includingthe ELODIE measurements (Fig. 8). The values of the goodness-of-fit measure ˜ l are notincreased very much with respect to those from Table 1 and are comparable with thosefrom Table 2. This means that one of the ACR fits can reflect the true configuration of thissystem quite well. The corresponding condition numbers C for the 2/1 resonance in Table 3are much less than in Table 1. This indicates that the topology of the likelihood surfacebecomes simpler and closer to the desirable paraboloidal one. The troubles connected withmultiple local maxima of the likelihood have been overcome in the ACR fits. We can nowsay more definitely, that the effect of the putative annual term in the Keck RV data on the -50 0-10050-50 0 50 100 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 R A D I A L VE L O C I T Y [ m / s ] -50 0-10050-50 0 50 100 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 R A D I A L VE L O C I T Y [ m / s ] Fig. 8
RV curves for theACR solutions I c A (toppanel) and I c B (bottompanel), plotted together withthe RV measurements andtheir residuals. The errorbars do not incorporateestimated values of theRV jitter. The best-fittingannual sinusoidal drift of theELODIE measurements waspreliminarily subtracted. fit quality is similar to the effect of the linear RV trend (which could be due to a long-period planet or brown dwarf in the system). These extra terms are significant in the fitsfor the 2/1 resonance (the ‘A’ group of solutions). For these fits, the estimations of falsealarm probabilities, calculated from the corresponding likelihood ratios according to Baluev(2008c), are about 2% and 0 .
3% (for the annual term and for the linear drift, respectively).The 5/2 resonance (B) solution does not require these terms in the RV model. This dilemmacan be solved by future observations only.
Now we consider the dynamical regimes of our nominal orbital configurations more closely.We are especially interested in how much the innermost planet can disturb the apsidal coro-tation of the two outer planets. This effect may be split in two categories:1. The planet ‘b’ can inspire some extra oscillation of the outer planets ‘c’, ‘d’ near theirunperturbed ACR (‘unperturbed’ means obtained without taking into account the influ-ence of the planet ‘b’).2. The planet ‘b’ can shift the position of the libration centre from the unperturbed ACR.Only the first effect may significantly affect the system stability, because only a large-amplitude oscillation around the libration centre can significantly increase the probabilityof close approaches of the resonant planets.Fig. 9 illustrates both effects. We can see that the orbital configuration I c A taken fromTable 3 shows moderate oscillation of the eccentricities e c , e d and libration of resonant an-gles s c , s d around some equilibrium values. However, the centres of oscillations are some-what shifted with respect to the unperturbed ACR. Heuristically, to counterbalance the grav-itational influence of the perturbing planet ‘b’, we need to increase somewhat the mass of E CC E N T R I C I T I ES o -60 o -30 o o o o
90 0 20 40 60 80 100 d s c s R ES O N A N T A N G L ES E CC E N T R I C I T I ES o -60 o -30 o o o o
90 0 20 40 60 80 100 d s c s R ES O N A N T A N G L ES E CC E N T R I C I T I ES o -90 o -60 o -30 o o o o o o
150 0 20 40 60 80 100 c s d sTIME [Kyr] R ES O N A N T A N G L ES Fig. 9
Temporal evolution of the ACR configurations of planets in HD37124 (for RV model I only). Leftpanels show the evolution of the eccentricities, right panels show the evolution of resonant angles. Top pair:resonance 2/1, orbital parameters were taken from Table 3, solution I c A. Middle pair: the same case, but thevalue of the mass m c was increased so that the value of ˜ K c increased by 1 m/s. Bottom pair: resonance 5/2(Table 3, solution I c B). The same character of motion is conserved on the time scale of 10 yr (and probablymuch longer) for all the cases. the outermost planet ‘c’. This assumption is confirmed by numerical integration: orbital con-figuration with m c increased by about 8% shows much less libration amplidute. Moreover,this adjustment decreases the statistic ˜ l by about 0 . c B), thelibration amplitude is quite small for the unperturbed ACR solution already. The evolutionof these orbital configurations appears perfectly regular and does not show any instability atthe time scale of at least 10 yr.More detailed analysis shows a great diversity of dynamical behaviour in the vicinityof the ACR solutions. The amplitude and character of the librations can be different. Forcertain orbital configurations inside the 2/1 MMR, the system may switch (from time totime) between alternating asymmetric ACRs with w c − w d ≈ − ◦ and w c − w d ≈ + ◦ .Librations surrounding simultaneously the pair of stable asymmetric stationary solutionsand unstable symmetric aligned ACR (see the left-middle panel in Fig. 7) are also possible. PE R I O D OG R A M VA L U E PE R I O D OG R A M VA L U E Fig. 10
The likelihood ratio periodograms of the RV residuals for different three-planet orbital models of thesystem of HD37124. Every value of these periodograms represents the modified likelihood ratio statistic (seeBaluev, 2008c), calculated for the base model of the residuals (free constant velocity offsets + free commonlinear drift) and for the alternative model incorporating also a sinusoidal variation with free amplitude andphase. The panels in the left column represent the periodograms constructed from the residuals of full RVdataset. They were constructed in the range of periods starting from 10 days, since the errors in the datesof ELODIE measurements do not allow to fit accurately more short periods. The graph in the right columnsshow the periodograms constructed for the Keck residuals only. The top pair of panels is for the ACR solutionin the 2/1 MMR (RV model I). The bottom pair of panels is for the ACR solution in the 5/2 MMR (also RVmodel I). For the case of graphs to the left, the normalized frequency bandwidth W ≈ W ≈ When the best fitting orbital structure of the planetary system appears unstable, we cansuspect that more planets orbit the star. It was the case for the system of m Ara: the de-termination of realistic orbits of m Ara b and c represented an essential difficulty (e.g.,Go´zdziewski et al., 2005) until the discovery of the planet m Ara e placed all in the places(Pepe et al., 2007; Go´zdziewski et al., 2007). To test the hypothesis of an extra planet or-biting HD37124, we use the likelihood ratio–based periodogram described in the paper(Baluev, 2008c). This periodogram represent a generalization of the usual Lomb (1976)-Scargle (1982) periodogram (as well as of the data-compensated discrete Fourier transformperiodogram by Ferraz-Mello (1981)) and incorporates a built-in estimation of the RV jitter.The graphs of such periodograms of RV residuals for several orbital solutions and in-volving different datasets are shown in Fig. 10. We can see that none of the peaks risesclearly beyond the apparent noise level. The formal significances of a periodogram peak canbe assessed using the analytic expression for the associated false alarm probabilities (i.e.,the probability to claim that the peak is statistically significant when actually it is a result ofnoise fluctuations) from the paper (Baluev, 2008a): ( false alarm probability ) ≈ We − z √ z , (13)where z is the height the maximum periodogram peak, and W = f max T eff is the normalizedfrequency bandwidth, with f max being the maximum frequency being scanned and T eff being the effective time series span (which is usually close to the actual time base). The applica-bility of this formula to the likelihood ratio periodograms was also discussed in (Baluev,2008c). None of the peaks on the periodograms in Fig. 10 possesses the false alarm proba-bility estimation less than ∼ c A (the upper-right panel in Fig. 10). This peak isclose to the period of 26 . ∼ P ≥
10 days (instead of P ≥ P d / ≈
450 days or P c / ≈
650 days would be extremelydifficult to extract from the synthetic RV curve. Such RV oscillation could be almost equallytreated as a Fourier overtone associated with Keplerian oscillations induced by the outerplanets (thus resulting in some change of their best-fitting eccentricities). Since the cur-rently available RV data for HD37124 do not allow reliable determination of the orbits inthe system even for the three-planet configuration, the more complicated (and worse deter-mined) four-planet configurations were not considered here. It is worth emphasizing thatthere is no necessity to call for four-planet configurations of HD37124. All present difficul-ties connected with this system can be explained from the positions of data analysis, that isby the lack of the RV data, which is not sufficient to constrain reliably the non-linear triple-planet RV model having large number of degrees of freedom. For example, it was shown byBaluev (2008b) that the unrealistically large formal estimations of the eccentricities can beinterpreted as a result of their statistical (systematic) biasing.
In the paper, the full set of high-precision RV data available for the planetary system ofHD37124 is analysed. The analysis involves different RV models and accounts for the re-quirement of the dynamical stability of the planetary system. The most likely orbital config-urations of the system, found in the paper, are split in four classes:1. Two outer planets ‘c’,‘d’ are captured in the 2/1 mean-motion resonance and move onorbits with low or moderate eccentricities ( e . . ever, these solutions show larger scatter of the RV residuals and thus are less likely thensolutions from other branches. When the eccentricies are fixed at small values, the ratioof the best fitting orbital periods of the two outer planets appears within a few per centof the 2/1 MMR. The values of the period ratio P c / P d exceeding 2 . − . P c / P d < . m Ara, which are close to, but likelynot trapped in, the 2/1 MMR with P b / P e ≈ . Acknowledgements
This work was supported by the Russian Foundation for Basic Research (Grant 06-02-16795) and by the President Grant NSh-1323.2008.2 for the state support of leading scientific schools. Iam grateful to Profs. K.V. Kholshevnikov and V.V. Orlov for useful comments and correction. I would liketo thank the referees, C. Beaug´e and the anonymous one, for careful reading of the manuscript and usefulsuggestions which helped to improve it.
References
Baluev R.V.: Assessing the statistical significance of periodogram peaks. Mon. Not. R. As-tron. Soc. , 1279-1285 (2008a)Baluev R.V.: Several problems of exoplanetary orbits determination from radial velocityobservations. In: Sun Y.-S., Ferraz-Mello S., Zhou J.-L. (eds.) Exoplanets: Detection,Formation and Dynamics (IAU Symp. 249). Suzhou, October 2007. Camb. Univ. Press,IAU Proc., , pp. 101-110 (2008b)Baluev R.V.: Accounting for velocity jitters in planet search surveys. Mon. Not. R. Astron.Soc., submitted, arXiv: 0712.3862 (2008c)Baluev R.V.: Optimal strategies of radial velocity observations in planet search surveys.Mon. Not. R. Astron. Soc., , 1375-1382 (2008d)Beaug´e C., Ferraz-Mello S., Michtchenko T.A.: Extrasolar planets in mean-motion reso-nance: apses alignment and asymmetric stationary solutions. Astrophys. J. , 1124-1133 (2003)Beaug´e C., Michtchenko T.A., Ferraz-Mello S.: Planetary migration and extrasolar planetsin the 2/1 mean-motion resonance. Mon. Not. R. Astron. Soc. , 1160-1170 (2006)Beaug´e C., Giuppone C.A., Ferraz-Mello S., Michtchenko T.A.: Reliability of orbital fits forresonant extrasolar planetary systems: the case of HD82943. Mon. Not. R. Astron. Soc. , 2151-2160 (2008)Butler R.P., Marcy G.W., Vogt S.S., Fischer D.A., Henry G.W., Laughlin G., Wright J.T.:Seven new Keck planets orbiting G and K dwarfs. Astrophys. J., , 455-466 (2003)Ferraz-Mello S.: Estimation of periods from unequally spaced observations. Astron. J. ,619-624 (1981) Ferraz-Mello S., Michtchenko T.A., Beaug´e C.: The orbits of extrasolar planets HD82943 cand b. Astrophys. J. , 473-481 (2005)Ford E.: Adaptive scheduling algorithms for planet searches. Astron. J. , 1008-1020(2008)Froeschl´e C., Lega E., Guzzo M.: Analysis of the chaotic behaviour of orbits diffusing alongthe Arnold web. Celest. Mech. Dyn. Astron. , 141-153 (2006)Gerassimov I.A., Vinnikov E.L., Mushailov B.R.: Canonical Equations in Celestial Mechan-ics [in Russian]. Moscow Univ. Press, Moscow (1996)Go´zdziewski K., Konacki M., Maciejewski A.J.: Orbital solutions to the HD160691( m Arae) Doppler signal. Astrophys. J. , 1136-1148 (2005)Go´zdziewski K., Konacki M., Maciejewski A.J.: Orbital configurations and dynamical sta-bility of multiplanet extrasolar systems around Sun-like stars HD202206, 14 Herculis,HD37124, and HD108874. Astrophys. J. , 688-703 (2006)Go´zdziewski K., Maciejewski A.J., Migaszewski C.: On the extrasolar multiplanet systemaround HD160691. Astrophys. J. , 546-558 (2007)Go´zdziewski K., Breiter S., Borczyk W.: The long-term stability of extrasolar system HD37124. Numerical study of resonance effects. Mon. Not. R. Astron. Soc. , 989-999(2008)Hadjidemetriou J.D.: Symmetric and asymmetric librations in extrasolar planetary systems:a global view. Celest. Mech. Dyn. Astron. , 225-244 (2006)Hadjidemetriou J.D.: On periodic orbits and resonance in extrasolar planetary systems. Ce-lest. Mech. Dyn. Astron., doi: 10.1007/s10569-008-9119-8 (2008)Ji J., Kinoshita H., Liu L., Li G., Nakai H.: The apsidal antialignment of the HD82943system. Celest. Mech. Dyn. Astron. , 113-120 (2003)Lee M.H., Peale S.J.: Secular evolution of hierarchical planetary systems. Astrophys. J. ,1201-1216 (2003)Lehman E.L.: Theory of Point Estimation. Wiley, New York (1983)Lomb N.R.: Least-squares frequency analysis of unequally spaced data. Astrophys. & Sp.Sci. , 447-462 (1976)Lissauer J.J., Rivera E.J.: Stability analysis of the planetary system orbiting u Andromedae.II. Simulations using new Lick Observatory fits. Astrophys. J. , 1141-1150 (2001)Michtchenko T.A., Beaug´e C., Ferraz-Mello S.: Stationary orbits in resonant extrasolar plan-etary systems. Celest. Mech. Dyn. Astron. , 411-432 (2006)Pepe F., Correia A.C.M., Mayor M., Tamuz O., Couetdic J., Benz W., Bertaux J.-L.,Bouchy F., Laskar J., Lovis C., Naef D., Queloz D., Santos N.C., Sivan J.-P., Sos-nowska D., Udry S.: The HARPS search for southern extra-solar planets. VIII. m Arae, asystem with four planets. Astron. & Astrophys. , 769-776 (2007)Scargle J.D.: Studies in astronomical time series analysis. II - Statistical aspects of spectralanalysis of unevenly spaced data. Astrophys. J. , 835-853 (1982)Short D., Windmiller G., Orosz J.A.: New solutions for the planetary dynamics inHD160691 using a Newtonian model and latest data. Mon. Not. R. Astron. Soc. ,L43-L46 (2008)Udry S., Mayor M., Queloz D.: Extrasolar planets: from individual detections to statisticalproperties. In: Deming D. and Seager S. (eds.) Scientific frontiers in research on extrasolarplanets. Washington, June 2002. ASP Conf. Series. vol. 294. pp. 17-27 (2003)Vogt S.S., Marcy G.W., Butler R.P., Apps K.: Six new planets from the Keck precisionvelocity survey. Astrophys. J. , 902-914 (2000)Vogt S.S., Butler R.P., Marcy G.W., Fischer D.A., Henry G.W., Laughlin G., Wright J.T.,Johnson J.A.: Five new multicomponent planetary systems. Astrophys. J. , 638-658 (2005)Voyatzis G., Hadjidemetriou J.D.: Symmetric and asymmetric 3:1 resonant periodic orbitswith an application to the 55Cnc extra-solar system. Celest. Mech. Dyn. Astron. , 259-271 (2006)Wright J.T.: Radial velocity jitter in stars from the California and Carnegie planet search atKeck observatory. Publ. Astron. Soc. Pacific , 657-664 (2005) A Obtaining ACR fits
Here we describe the procedure of obtaining the ACR fits in more details. The algorithm that we are aboutto describe may be useful also for other planetary systems. But before we proceed, we have to choose somecoordinate system. Actually, in the case of HD37124, it may be checked that the offsets of resulting orbitalparameters, referenced in different coordinate systems, would be quite negligible in the sense of the RV fitquality. However, the type of the coordinate system should be stated for the purposes of long-term integrationsof the planetary system: the coordinate system used in the integration should match the given orbital elements.The actual choice of the Jacobi coordinates was motivated here by the fact noted in (Lissauer & Rivera, 2001;Lee & Peale, 2003), that it is Jacobi coordinate system in which the osculating orbital elements are mostlyclose to those obtained using the kinematic (Keplerian) RV model, especially when the system containshierarchical planet pairs (like pairs of planets b–c and b–d in the case of HD37124). This property of theJacobi coordinates may be useful, for instance, during a transition from multi-Keplerian to N-body modelof the RV, that was (and probably will be) needed for some resonant planetary systems after accumulating asufficiently long observation time span.In the Jacobi coordinates, the Hamiltonian of a system with N planets looks like H = N (cid:229) i = p ′ i m ′ i − Gm m i r i ! − (cid:229) ≤ i < j ≤ N Gm i m j r ij , (14)where G is the gravitational constant, m i are the planetary masses, r i are the astrocentric distances of theplanets, r ij are the distances between the planets, p ′ i are the Jacobi momenta and m ′ i = m i ( (cid:229) i − j = m j ) / ( (cid:229) ij = m j ) are the Jacobi masses. Now we need to split (14) in the Keplerian part that we consider as unperturbed oneand in the part that we consider as perturbational function. This may be done non-uniquely. We adopt thesplitting H = (cid:229) Ni = H Kep , i − (cid:229) Ni = R i , where H Kep , i = p ′ i m ′ i − GM i − r ′ i m i ! , R i = Gm i N (cid:229) j = i + m j r ij + m r i − M i − r ′ i ! . (15)Here r ′ i is the length of the Jacobi radius-vector for the i th planet, M i − = (cid:229) i − j = m j is the sum of the starmass and of the masses of all planets being interior with respect to the given planet. Therefore, H Kep , i ischosen so that the corresponding unperturbed Keplerian orbit is referenced to a fictitious central body havingmass M i − . We will use only the second-order approximation of the Hamiltonian. In this approximation, werepresent R i = (cid:229) Nj = i + R ij + O ( m ) , where R ij = Gm i m j | r ′ i − r ′ j | − r ′ i · r ′ j r ′ j − r ′ j ! . (16)Here r ′ i are the Jacobi position vectors. This approximation for the case of two planets may be found in(Gerassimov et al., 1996, § N ≥ R ( a , e inn , e out , l inn , l out , Dw ) = | r inn − r out | − r inn · r out r − r out . (17)Here, the vectors r inn and r out describe positions of abstract ‘inner’ and ‘outer’ planets moving on Keplerianorbits. The Keplerian orbital elements marked by the subscript ‘inn’ refer to the orbital elements of the innerplanet, and quantitties marked by the subscript ‘out’ correspond to the outer planet. The semi-major axis of7the inner abstract planet is set to unit, and the semi-major axis of the abstract outer planet is equal to a . Theparameter Dw = w out − w inn represents the difference between the longitudes of periastra. Then the obviousidentity R ij = Gm i m j a i R ( a j / a i , e i , e j , l i , l j , w j − w i ) holds true.As usually, the Keplerian parts of the Hamiltonian depend only on the Delaunay actions L i (and, as onparameters, on the masses of the planets) and remain unchanged by any averaging. Let us now restrict ourattention to the motion of the (possibly) resonant planets HD37124 c and d, neglecting the influence of theinnermost planet b. The perturbational function describing the interaction of the planets c and d is givenby Gm c m d a d R ( a = a c / a d , e d , e c , l d , l c , Dw = w c − w d ) . The corresponding averaged perturbational functionis then given by Gm c m d a d h R i ( a , e d , e c , s d , s c ) , where the averaged function h R i depends on only five inputparameters a , e inn , e out , s inn , s out : h R i ( a , e inn , e out , s inn , s out ) = Z p R ( a , e inn , e out , s inn + p q , s out + q q , s inn − s out ) d q p . (18)Note that the last term in the expression (17) is averaged to 1 / a . This term reflects only the fact that theosculating Keplerian orbits are refered to different fictitious central masses. In fact, we need to average onlythe resting classical expression of the perturbational function: h R i = (cid:28) | r inn − r out | − r inn · r out r (cid:29) − a . (19)We perform the averaging (18,19) by means of numerical integration tools, as it was proposed by Michtchenko et al.(2006). This way of numerical averaging of the Hamiltonian is very easy to implement and simultaneously isquite rapid and precise. In the same way, we can calculate various derivatives of the averaged function h R i ,that we will need below. For example, a ¶ h R i ¶a = (cid:28) r inn · r out − r | r inn − r out | + r inn · r out r (cid:29) + a , (20)where the averaging can be again performed numerically.In the next step, we need to build in our RV fitting algorithm the four equality bounds ¶ h H cd i / ¶ ( I c , I d ) = ¶ h H cd i / ¶ ( s c , s d ) = ¶ h R cd i / ¶ ( s c , s d ) =
0, which determine the location of the ACR. The secondpair of equations reflects the fact that s d and s c should correspond to the extremum value of h R i ( a = a c / a d , e d , e c , s d , s c ) (considering the values of e cd and a fixed). This condition can be used to construct theangles s c , d as functions of the eccentricities e c , d : s d = s ∗ inn ( a , e d , e c ) and s c = s ∗ out ( a , e d , e c ) , where the func-tions s ∗ inn ( a , e inn , e out ) and s ∗ out ( a , e inn , e out ) provide an extremum to h R i ( a , e inn , e out , s inn , s out ) given fixed a , e inn , e out . When the ACR is symmetric, the functions s ∗ inn and s ∗ out can be found easily (actually, they ap-pear to be constant in this case). However, for an asymmetric ACR, we have to locate the values of s c and s d numerically. Eventually, this numerical procedure is still sufficiently rapid and can be implemented as a‘black-box’ subroutine, which returns the values of s ∗ inn and s ∗ out for input values of e inn , e out , and a . Note thatsince we consider here only MMR solutions, a ≈ a ≡ ( p / q ) / with a error of O ( m c , d / M ⋆ ) . When calcu-lating the resonant angles, we can quite neglect such errors and put simply a = a . Such errors in resonantangles will not produce significant changes in the dynamics of the planetary system. The resulting values of s c and s d and the definitions (10) can be used in the work of the fitting algorithm to express w c and w d viathe other free parameters: l c , l d , e c , e d .Another pair of constraints, ¶ h H cd i / ¶ ( I c , I d ) =
0, can be transformed to a more convenient (and equiv-alent) form, involving partial derivatives over the Delaunay actions L c , d and G c , d = L c , d h c , d (here h c , d = − e c , d ). The first transformed equation, ¶ h H cd i / ¶ G c = ¶ h H cd i / ¶ G d , reflects the coincidence of the seculardrifts of the orbital periastra. After neglecting insignificant errors of the order of the planetary masses, thisequation can be simplified to˜ K d ˜ K c = a r ( a , e d , e c ) with r ( a , e inn , e out ) = ¶ h R i / ¶h inn ¶ h R i / ¶h out (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s inn = s ∗ inn ( e inn , e out ) s out = s ∗ out ( e inn , e out ) , (21)where we note that ¶¶h = − h e ¶¶ e . The equality (21) can be used to express the RV semi-amplitude of one ofthe resonat planets via the RV semi-amplitude of another one and via the orbital eccentricities. The secondtransformed equation, p ¶ h H cd i / ¶ L c = q ¶ h H cd i / ¶ L d , reflects the vanishing of the secular drift of the critical8angle pl c − ql d (where l are the planetary mean anomalies). This should provide the long-term constancy ofthe planetary conjunction positions. After some simplifications, this equation may be rewritten as P d P c pq = − m c M ⋆ n ( a , e d , e c ) with n ( a , e inn , e out ) = (cid:20)(cid:18) pq h out − h inn (cid:19) ¶ h R i ¶h inn − (cid:18) a ¶ h R i ¶a (cid:18) + pq r (cid:19) + h R i (cid:19)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s inn = s ∗ inn ( e inn , e out ) s out = s ∗ out ( e inn , e out ) . (22)This equation introduces a error of the second order only, that is O (cid:0) ( m / M ⋆ ) (cid:1) . Actually, specifically to thesystem of HD37124, the small O ( m / M ⋆ ) deviation of the period ratio from the exact resonance does not affectsignificantly the quality of the RV fit (moreover, the values of n for the ACR fits from Table 3 appeared lessthan 0 . P c = pq P d .However, it is this O ( m / M ⋆ ) period ratio deviation that determines the secular drift of the critical angle pl c − ql d and the ACR state of the system. Therefore, in general case it is necessary to take this deviationinto account when constructing the ACR fits. The equation (22) can be used to express the osculating orbitalperiod P c or P d via the resting free variables. Note that the deviation of the ratio of the osculating semi-majoraxes can be written down (to within the first order) as aa = + m c M ⋆ n ( a , e d , e c ) + . (23)Therefore, for any given values of the parameters e c , e d , l c , l d , ˜ K d , P d we can obtain the ACR val-ues of the parameters w c , w d , K c with a error of O ( m c , d / M ⋆ ) and the ACR value of P c with a error of O (cid:0) ( m c , d / M ⋆ ) (cid:1)(cid:1)