Orbital structure of the GJ876 extrasolar planetary system, based on the latest Keck and HARPS radial velocity data

# Orbital structure of the GJ876 extrasolar planetary system, based on the latest Keck and HARPS radial velocity data

Roman V. Baluev Pulkovo Astronomical Observatory, Pulkovskoje sh. 65/1, Saint Petersburg 196140, Russia;
Sobolev Astronomical Institute, Saint Petersburg State University, Universitetskij pr. 28, Petrodvorets, Saint Petersburg 198504, Russia
11email: roman@astro.spbu.ru
###### Abstract

We use full available array of radial velocity data, including recently published HARPS and Keck observatory sets, to characterize the orbital configuration of the planetary system orbiting GJ876. First, we propose and describe in detail a fast method to fit perturbed orbital configuration, based on the integration of the sensitivity equations inferred by the equations of the original -body problem. Further, we find that it is unsatisfactory to treat the available radial velocity data for GJ876 in the traditional white noise model, because the actual noise appears autocorrelated (and demonstrates non-white frequency spectrum). The time scale of this correlation is about a few days, and the contribution of the correlated noise is about  m/s (i.e., similar to the level of internal errors in the Keck data). We propose a variation of the maximum-likelihood algorithm to estimate the orbital configuration of the system, taking into account the red noise effects. We show, in particular, that the non-zero orbital eccentricity of the innermost planet d, obtained in previous studies, is likely a result of misinterpreted red noise in the data. In addition to offsets in some orbital parameters, the red noise also makes the fit uncertainties systematically underestimated (while they are treated in the traditional white noise model). Also, we show that the orbital eccentricity of the outermost planet is actually ill-determined, although bounded by . Finally, we investigate possible orbital non-coplanarity of the system, and limit the mutual inclination between the planets b and c orbits by , depending on the angular position of the mutual orbital nodes.

###### Keywords:
extrasolar planets radial velocity red noise mean-motion resonance
journal: Celestial Mechanics & Dynamical Astronomy

## 1 Introduction

The first extrasolar planet in the system orbiting the red dwarf GJ876 was detected independently by Delfosse et al. (1998) and Marcy et al. (1998) on the basis of high-precision radial velocity (hereafter RV) Doppler measurements, acquired at ELODIE+CORALIE and Keck HIRES spectrographs. Soon after this, the second planetary companion was discovered by Marcy et al. (2001). That discoveries showed that the exoplanetary system of GJ876 is an extraordinary object. The two planets b and c were massive giants (having minimum masses of roughly and of Jupiter masses) orbiting in short period orbits (approximately and days). Therefore, it became clear that the planets orbital dynamics should be significantly affected by the 2/1 mean-motion resonance (hereafter MMR). The resonance itself is not yet very astonishing, because now we already know many extrasolar planets in various MMRs. The extraordinarity of the GJ876 system comes from the fact that interplanetary gravitational perturbations were directly detected in the observed radial velocity data, i.e. they reveal themselves on the observational time scale (Laughlin & Chambers, 2001; Rivera & Lissauer, 2001).

Still, no other extrasolar planetary system around a main sequence star is known to demonstrate so clear and well-measurable signatures of planetary perturbations in its radial velocity time series (although, there are some exoplanetary systems where such perturbations are suspected to be non-negligible). A few favoring factors meet each other in GJ876: large masses of the planets b and c, small star mass (hence especially large planet/star mass ratios), rather short orbital periods (thus shorter perturbations time scale), and, of course, the low-order MMR. The dynamical perturbations reveal themselves mainly in the form of secular circulation of the planetary apsidal lines, triggering slow change in the non-sinusoidal shape of the RV oscillations being observed. The notion “secular”, however, looks somewhat odd here, since the period of this “secular” circulation is only about  yrs. During the whole observation term since 1998 till today, these apsidal lines completed roughly a single revolution each.

This dynamical effect complicates the procedure of the RV data analysis, but in exchange it enables us to determine (solely from the RV data) the inclination of the system to the sky tangent plane. This allows us to determine, instead of the minimum planetary masses , the true masses , which normally remain unconstrained in the RV exoplanet detections.

Later, Rivera et al. (2005) reported the discovery of the third planet in the system, based on further Keck RV observations of GJ876. This planet (d) possesses very low mass of Earth masses and a very short orbital period of approximately  days. In the paper (Bean & Seifahrt, 2009), the question of orbital non-coplanarity of this system was studied extensively, and these authors gave an estimation of for the mutual orbital inclination of the planets b and c. This result was based on the Keck RV data from (Rivera et al., 2005) and on the HST astrometry data from (Benedict et al., 2002). Recently, Correia et al. (2010) reanalyzed these old Keck data adding to them a new array of very accurate RV measurements, obtained by the HARPS spectrograph (ESO). He presented improved three-planet orbital fits, which generally agree with the fits by Rivera et al. (2005).

Finally, in the very recent paper (Rivera et al., 2010) the discovery of the fourth Uranus-mass planet e was reported. This last planet is also remarkable, because it forms a Laplace three-planet resonance with two other giant planets, so that the ratio of the periods is close to .

It is interesting, however, that the RV signature of this fourth planet had been found in the old Keck data already, e.g. in (Baluev, 2008b). It is notable that the fourth planet parameters from (Baluev, 2008b) are surprisingly close to the ones given by Rivera et al. (2010), including even the mean longitude and the small orbital eccentricity. This orbital configuration appears pretty stable at long time scales. Rivera et al. (2010) also mention that they were seeing the signs of the fourth planet in their data since 2005, but until 2010 they could not find a stable four-planet solution, due to the large planet e eccentricity. The reason of such difference probably comes from the fact that Rivera et al. (2010) did not take into account the annual errors of a few m/s, which were detected in the old Keck data by Baluev (2008b). Such annual errors may appear relatively frequently in the RV planet search surveys, as discussed in (Baluev, 2009). They may have different physical sources; in case of GJ876 they were probably caused by some errors in the data reduction pipeline, since new Keck data from Rivera et al. (2010) look free from such errors. In fact, a careful error analysis of the RV data could enable the robust detection of the fourth planet several years ago already.

Since the quality of the Keck data is considerably improved in (Rivera et al., 2010), and the completely new very accurate RV data is now published (Correia et al., 2010), we conclude that it is time to perform such detailed analysis now. Accurate and unbiased values of planetary masses and orbital parameters in so remarkable system as GJ876 represent a huge interest for the celestial mechanics studies, as well as for more astrophysical branches like the tidal star-planet interaction studies (note the “hot super-earth” GJ876 d) or for the planet formation theories. In particular, it is important to have some reliable estimations or limits of the orbital non-coplanarity of such system, or to know how large the orbital eccentricity of the “hot Earth” GJ876 d actually can be.

The plan of the paper is as follows. First, in Sect. 2 we describe in detail the fast approach of fitting the dynamically perturbed planetary configurations on the basis of the star RV time series. We also explain the conventions we adopt to determine the reference osculating orbital parameters, and describe the numerical integrator we use in this study. In Sect. 3 we describe in more details the RV data we use in the paper and also give some preliminary analysis results, obtained during direct fitting of these data. In Sect. 4, we find that these RV data cannot be modeled in a traditional way, because they contain significant fraction of autocorrelated (non-white) noise. Also, we propose here a modified RV fitting algorithm, which allows to take such correlated RV errors into account properly. In Sect. 5 we study the planet e orbital parameters and show that its eccentricity is still ill-determined although bounded by from the upper side. In Sect. 6 we give final estimations of all planetary parameters, taking into account the effects of the RV noise correlateness and of the bad determinability of the planet e eccentricity. In Sect. 7, we study the (non-)coplanarity of the system in view of currently available RV data. Finally, Sect. 8 describes briefly the dynamical evolution of the admissible GJ876 orbital configurations.

## 2 Fitting graviationally perturbed orbits

### 2.1 Evaluating radial velocity function

Before we proceed to the investigation of the GJ876 planetary system itself, we describe the algorithm we used to perform -body fitting of RV data. Although such -body fitting have been already carried out in several previous works (Laughlin & Chambers, 2001; Rivera & Lissauer, 2001), but they omit important details of the algorithms used. We may highlihgt a very recent paper (Pál, 2010), which introduces a method of orbital fitting, which is similar (in spirit) to the method that we use here. However, many ideas of the methods developed by Pál (2010) look pretty different, so we still need to describe our method in detail here. The algorithm that we use represents a variation of the popular non-linear Levenberg-Marquardt minimization algorithm with a minor modification necessary to take into account the “RV jitter” phenomenon, as discussed in (Baluev, 2009) and below. This algorithm is an iterative gradient method, which requires to evaluate (on different iteration stages) the fittable model (RV curve model, in our case) and its partial derivatives over all free parameters. In case of GJ876, the two latter subtasks are non-trivial, since we need to take into account the full dynamical model of planetary perturbations.

Let us write down the equations of the -body problem in the astrocentric coordinate system:

 1k2d2→ridt2=−(1+μi)→rir3i+∑j=1..Nj≠iμj(→rj−→ri|→rj−→ri|3−→rjr3j),i=1,2,…,N. (1)

Here are the ratios of the planet masses to the star mass , and is the analogue of the Gauss gravitational constant. Note that further we may call just “planetary masses” for shortness, since usually it should not introduce any misunderstandings. We assume that axis is directed along the observer’s line of sight, and two other axes are directed arbitrarily in the tangent sky plane.111Note that the term “-body”, which we use in the paper frequently, does not imply that the number of the bodies in the system is denoted by . Instead, stands for the number of interacting planets, for the number of RV measurements, and “-body” must be understood as a solid term.

The initial conditions for the system (1) may be expressed via some set of osculating orbital elements forming vectors using the well-known functions of the Keplerian motion. However, there is no unique definition of the osculating Keplerian elements. Their set depends not only on our choice of the parametrization, but also on the non-unique splitting of the forces in (1) into the unperturbed Keplerian and perturbational parts. For instance, we can assume that the Keplerian part of the force is , and everything else in the right hand side of (1) is the perturbation. In such a case, the initial planetary coordinates and velocities could be expressed via corresponding osculating orbital elements as

 →ri|t=0=k23→Kr(→πi),→vi|t=0=k23→Kv(→πi), (2)

where the functions and should represent the solution of the standard equation of the Keplerian motion at . Each vector incorporates six osculating orbital elements of the planets: orbital period , eccentricity , pericenter argument , mean argument of latitude (mean anomaly + ), inclination to the sky plane , longitude of the ascending node . The specific of the exoplanetary orbit determination problem sets a few pecularities of these parameters that we need to highlight:

1. We remind the reader that the angle is traditionally counted from the descending orbital node (Ferraz-Mello et al., 2005, sect. 2.2). Historically, this deviation from the classical definition of occured because exoplanet researchers traditionally model the host star’s RV wobble using the Keplerian formula , which actually expresses the RV of the orbiting planet itself. This shifts the derived value of by , due to the RV sign change. Alternatively, we may think that refers to the star’s orbit around the common barycenter, although such treatment becomes a bit vague for multi-planet systems. This pecularity of the angle is insignificant for the most of the exoplanets, since the positions of the planetary orbital nodes usually remain unconstrainable anyway. It may become important for GJ876, however.

2. Consequently, the angle is also counted from the descending node. We use instead of the full mean longitude , since the absolute positions of the orbital nodes still remain unconstrained until the astrometric observations are used. Radial velocities can only constrain the differences between . For coplanar configurations, dealing with and is equivalent to dealing with and the pericenter longitude , since all are equal to the same constant value.

3. In the equations (2), the dependence on the actual (in fact, on ) is extracted separately in the external factors . Such dependence on appears when contains planetary orbital period (or e.g. mean motion) as a primary parameter, and the semi-major axis is treated as only a derived one. We use such “period-oriented” choice because it is the planetary orbital period (or the period of corresponding RV oscillation) which is drawn from the observations directly, rather than the orbital semi-major axis. If we chose the semi-major axis as a primary parameter, the factors in (2) would be and .

The astrocentric coordinate system makes the -body equations (1) simple for numerical integration, but it is not well suitable to reference orbital elements. This is mainly because of relatively significant dependence of the planetary orbital periods on the choice of the coordinate system, which has been already discussed in (Lissauer & Rivera, 2001; Lee & Peale, 2003; Baluev, 2008c). As it was noted in these works, the osculating orbital periods (or mean motions) referenced in the Jacobi coordinate system generally match better the apparent orbital periods (or apparent average drift of the mean longitude). In practice that would mean that, e.g., we may freely experiment with turning on/off the gravitational interactions of any planet without the need to manually adjust its orbital period when switching from one model to another. Otherwise, we should remember that there may be a significant offset between the best fitting orbital period in the Keplerian RV model framework (which implies apparent orbital period) and the full -body framework (the osculating period). Although such offsets are small (), they may be comparable to or even significantly exceed the relative statistical error of the respective periods (e.g. in case of the innermost planet of GJ876). In such a case, the orbital fitting is slowed down and eventually may even fail to converge to the correct period value (it may be attracted by a spurious close alias or just noisy neighboring period). The choice of the Jacobi system usually eliminates these troubles, if we also assume that each osculating Keplerian orbit refers to a fictitious central mass incorporating the star’s and this planet’s mass and masses of all inferior planets.

According to (Baluev, 2008c, appendix A), one way to define the unperturbed part of the -body Hamiltonian in the Jacobi coordinates is

 HKep,i=γiγi−1→p′i22μi−γi−1k2μir′i,γi=1+i∑j=1μi. (3)

We adopt this (non-unique) way to split the Hamiltonian into the Keplerian and perturbational parts, because the corresponding equation of the Keplerian motion

 d2→r′idt2=−k′i2→r′ir′i3,k′i2=k2γi=GM⋆γi (4)

involves the value of the central mass instead of . This is exactly what we need. The initial conditions for the Jacobi coordinates and velocities should represent a solution of the latter equation and therefore should look like

 →r′i∣∣t=0=k′i23→Kr(→πi),→v′i∣∣t=0=k′i23→Kv(→πi) (5)

These formulae differ from (2) only by the coefficient , which now involves planetary masses too. Note that in case of a highly hierarchical system with negligible mutual perturbations between planets, such definition of the osculating orbital periods would make them infinitesimally close to the apparent revolution periods. The exoplanetary systems are not always hierarchical, but in practice (in particular, for GJ876) the offsets between apparent and so-defined osculating planet periods in usually remain satisfactory.

After that, we can transform the Jacobi vectors (5) to the astrocentric ones using the formulae

 →ri = →Ti(→r′1,…,→r′i,μ1,…,μi−1)≡→r′i+i−1∑j=1μjγj→r′j, →vi = →Ti(→v′1,…,→v′i,μ1,…,μi−1). (6)

The formulae (5) and (6) jointly express the Cartesian initial conditions in the astrocentric system via the osculating orbital elements defined in the Jacobi system, and via the planet masses . For shortness, we can write down these compound dependences as

 →ri|t=0 = →Ri(k,→π1,…,→πi,μ1,…,μi), →vi|t=0 = →Vi(k,→π1,…,→πi,μ1,…,μi). (7)

Note that these functions also involve the parameter , which depends on the star mass , which we assume is a priori known from astrophysical models of its spectrum. Following Correia et al. (2010), we assume throughout the paper.

Having the astrocentric initial conditions from (7), we can integrate (1) to the desired time. After that, we will have astrocentric planetary positions and velocities. We prefer to integrate the astrocentric motion equations (1), since the analogous equations in the Jacobi system would be significantly more complicated and thus would slow the algorithm down. After the integration, the barycentric velocity vector of the star can be obviously expressed via the planetary astrocentric velocities as

 →v⋆=−N∑i=1μi→vi/(1+N∑i=1μi). (8)

Finally, all stages explained in this section allow us to determine the model RV curve, based on the planetary orbital parameters and masses as fittable free variables. The total observable radial velocity of the star incorporates also a constant radial velocity (we denote it as ), caused by the motion of the planetary system barycenter. Moreover, since the RV measurements available for GJ876 are basically relative, we should also introduce different fittable values of for the data series obtained at different instruments.

### 2.2 Using variational equations

Probably the most important stage of the Levenberg-Marquardt algorithm (as well as of any other gradient optimization method) is evaluation of the partial derivatives of the data model over the free parameters. In problems with complicated model function (like -body RV fitting) most of the computational time is spent during evaluation of these derivatives. The simple way to evaluate them is just to use numerical differentiation of the original RV model (which should be calculated using -body integration). This way is easy to implement algorithmically, but this leads to a very bad speed/error ratio, since the error of numerical differentiation is considerably larger than the error of the original function. A better way to evaluate these derivatives is to integrate variational equations of the -body problem (in the literature on non-linear optimization they are also called sensitivity equations, see e.g. chapter 8 in the book by Bard (1974)). It allows to obtain roughly the same accuracy of derivatives as in the original function with the same integration step and number of equations to integrate. Planetary positional vectors and their velocities depend on the time, on the planetary masses and on the orbital parameters . Differentiation of (1) over and yields the following linear ODE system for partial derivatives of :

 1k2d2dt2∂→ri∂→πj = −(1+μi)Π(→ri)∂→ri∂→πj+ +∑k=1..Nk≠iμj[Π(→rj−→ri)(∂→rj∂→πj−∂→ri∂→πj)−Π(→rj)∂→rj∂→πj], 1k2d2dt2∂→ri∂μj = −→rjr3j+(1−δij)→rj−→ri|→rj−→ri|3−(1+μi)Π(→ri)∂→ri∂μj+ (9) +∑k=1..Nk≠iμj[Π(→rj−→ri)(∂→rj∂μj−∂→ri∂μj)−Π(→rj)∂→rj∂μj], Π(→a)=a2I−3→a⊗→aa5,δij={1,i=j0,i≠j.

Here is the Kronecker delta, is the identity matrix, denotes the dyadic product of vectors, and thus is a matrix having elements . The derivatives should be treated as Jacobian matrices. The system (9) should be integrated simultaneously with the original equations (1). The partial derivatives of are equal to the temporal derivatives of the corresponding partial derivatives of and can be obtained during the integration automatically. The initial values for and can be obtained by means of differentiation of the compound functions (7) over the corresponding parameters. We do not give here the resulting expressions since they are rather complicated and numerous, and actually their main parts involving derivatives of over can be found in the classical literature on orbit refinement (e.g. Duboshin et al., 1976, sect. 3.3). The derivatives of (7) over are also clearly calculatable, but they are rather complicated as well, and we have to omit them too.

### 2.3 Numerical integrator

We need some numerical integration method to solve the differential systems (1) and (9). We do not need it to be well-behaved (in any sense) on long time scales, since our integration timescale is rather short ( yrs or  revolutions of the inner giant planet). But we do need it to be fast on this timescale. After a few experiments, we suggest that the integrators of the Everhart type (Everhart, 1973, 1974) would meet our requirements. We say “Everhart type” since we actually used several integrators belonging to a general family of integrators similar to those constructed originally by Everhart. The original Everhart integrator was based on the Radau quadrature formula, based on certain asymmetric splitting of the integration step. In accordance with Avdyushev (2010)222See preprint at http://www.scharmn.narod.ru/AVD/Gauss_15_2.pdf, in Russian, it is possible to construct an integrator based on an arbitrary splitting of the integration step. These integrators represent, in fact, implicit Runge-Kutta integrators, equipped by an efficient method of evaluating the predictors on each step. Certain segment splittings (like the Radau one) allow to increase the order of the integrator significantly. Moreover, it is known that the use of the Legendre splitting makes the integrator symplectic. Symplectic integrators are very well suitable for long-term integrations, because they can preserve the Hamiltonian properties of the original -body problem over a much longer term (e.g., they show much slower energy error accumulation). It is important that in the case of the Everhart integrators the symplecticity can be obtained even with no significant trade-off or undesired side effects.

The only significant disadvantage of such integrators is that they require constant time step for symplecticity. The algorithm itself does contain an optional step-size control mechanism, but variable step destroys the symplectic property. Therefore, we choose to use the same symplectic integrator both for the short and long integration terms. This is the -th order Everhart type integrator, based on the -node Legendre splitting. For short-term integrations (needed to perform orbital fits) we use the variable step. For long-term integrations (stability tests), we hold the step size fixed.

The algorithm that we use is largely based on the Avdyushev’s FORTRAN program code with several modifications, aimed to increase its performance. We omit the detailed description here, since we plan to provide it in a separate work in the future, along with the program source code.

## 3 Radial velocity data for GJ876 and its preliminary planetary orbital solution

The main datasets available are Keck RV measurements from (Rivera et al., 2010), and HARPS RV measurements from (Correia et al., 2010). Correia et al. (2010) also use RV datasets obtained at ELODIE and CORALIE spectrographs, and there is also Lick RV dataset in (Marcy et al., 2001). We include these data too, although it is necessary to note that they are considerably less accurate and in fact have insignificant effect on the system orbital solution — only a few per cent of the parameters uncertainties. The average internal RV precision of the HARPS and Keck datasets is about and  m/s, their time span about and years, respectively. The average RV uncertainties of the Lick, ELODIE, and CORALIE datasets are about  m/s or more, and they span and years and months, respectively.

The internal errors stated in the published RV tables do not yet constitute the full RV uncertainty. It is well-known that high-precision RV exoplanet searches suffer from the phenomenon of the so-called RV “jitter”, which adds extra irregular variations in their RV data. This excessive jitter was initially explained via star astrophysical activity effects, but later it was shown that various systematic instrumental effects (Baluev, 2009) may be significant as well (and may even dominate sometimes). Anyway, this jitter has to be taken into account during any further analysis. It is important here that in practice the effective value of the apparent RV jitter of the same star is usually different for different instruments, implying different (and poorly known) weighting coefficients for different datasets. To take this jitter into account properly, we utilize the maximum-likelihood approach suggested in (Baluev, 2009). This approach is based on maximizing the datasets’ joint likelihood function, which depends on the RV curve (planetary) parameters as well as on the parameters determining the statistical structure of the errors in the RV data. This allows to estimate all necessary jitter values simultaneously with planetary parameters.

According to (Baluev, 2008c), we will obtain the best fitting values of all involved parameters by means of maximizing the objective function

 ln~L=−J∑j=1Nj∑i=1(lnσfull,ji+12γ(rjiσfull,ji)2)−Nln√2π, (10)

where stands for the number of the data points in the th dataset, denotes the number of the datasets ( in our case), is the total number of observations, is the th RV residual (O-C) in the th dataset, and is the full variance of the corresponding RV measurement, which incorporates the internal measurement variance as well as the jitter variance . The quantity (not to be mixed with gammas in (3), which are indexed) is equal to , where is the number of free parameters of the RV curve (number of degrees of freedom in the RV model). This correcting divisor helps to remove the systematic bias in the RV jitter estimations, as discussed in (Baluev, 2009). This is necessary, because the residuals to the best fitting model are always systematically smaller than the original data errors, and without any correction they would yield underestimated value for the jitter.

The function (10) depends on the RV curve parameters (via , which involve the RV model) and on the RV jitter values (via ). The values of these parameters, where the function (10) reaches its maximum, represent the necessary best fitting estimations. To estimate a quality of a given fit, we use the statistic defined in (Baluev, 2009). This quantity represents a monotonous function of , but is measured in m/s and thus is intuitively more clear and comparable to more traditional measures like r.m.s. (which we will use too, for some comparison).

Other details of this method can be found in (Baluev, 2008c, 2009). Here we use the same formalism, with the major difference that we now deal with the full gravitational -body RV model, instead of the multi-Keplerian one. Also, on contrary with (Baluev, 2008c), we do not add in the RV model the terms describing possible annual systematic errors in the data. Although such annual errors existed in the old Keck data from (Rivera et al., 2005), it seems the new data should have been fully corrected (see Fig. 1). The HARPS observations for other stars seemingly were always free from such errors.

For each planet, we will fit the following osculating orbital parameters: orbital period (equiv. to the mean motion), mean argument of latitude , eccentricity , pericenter argument . We also fit common orbital inclination (assuming coplanar system). In addition to the orbital parameters, we should also fit the planet masses. However, it is not convenient to adopt planet masses as primary fit parameters. In the traditional non-perturbed framework, the masses themselves are not determinable at all. In such a case, the RV oscillation semi-amplitude, , is fitted (as primary parameter), which allows to derive the minimum planet mass . It is well-known that without detectable interplanetary gravitational perturbations the inclination remains unconstrained, and the true planet mass remains unknown. In our case the system inclination is well-constrainable, but only for coplanar model. We will consider non-coplanar configurations below too, where individual inclinations may be poorly determined. In such a case it would be better not to mix the uncertainties of the RV amplitude and inclination, otherwise we risk to deal with significant troubles during numerical fitting, due to strong correlations between various fit parameters. Also, we usually will not include the innermost planet in the -body integration, assuming that on the observation timescale its motion is close to a Keplerian orbit. For this planet, we just have to leave with the RV semi-amplitude.

Therefore, it would be better still to adopt the RV semi-amplitudes as primary fit parameters, instead of the planet masses. But then we must clarify what is “RV semi-amplitude” in the perturbed case. Actually, we cannot strictly define any “amplitude” for the perturbed motion, since this motion is not periodic. Nevertheless, in the case of GJ876 the deviations from the strict periodicity are small, and we still can determine the RV amplitude approximately, within the error of . Eventually, we need just to bind this RV amplitude parameter to the corresponding planet mass. Therefore, we can just define it via using some simple formula close in shape to the formula from the Keplerian RV case. We adopt the following definition:

 μ=K√1−e2sini(P2πk2)1/3. (11)

We must emphasize that this formula no longer expresses planet mass via the corresponding RV semi-amplitude, as it would be in the case of unperturbed motion (). Rather, it now defines the “RV semi-amplitude” via the planet mass. We may think of (11) as of a definition of the osculating RV semi-amplitude, which completes the set of usual osculating orbital elements. So-defined “osculating RV semi-amplitude” is in fact just an intermediate fit parameter needed to separate the uncertainty of the orbital inclination from the uncertainty of the planet mass. The apparent semi-amplitude of the star’s RV oscillation, caused by the corresponding planet, should be very close to the value of defined in (11). Furthermore, we decide to use not even the semi-amplitude itself, but the quantity , as it was done in (Baluev, 2008c). Such choice allows to eliminate the eccentricity from the relation (11) and thus facilitates the conversions between and .

To obtain some basic preliminary estimations of the system parameters, we take the Keck-only orbital solution from (Rivera et al., 2010) as a starting approximation and perform the non-linear maximization of the likelihood function, as it was explained above. Table 1 contains the resulting best fitting estimations of all planetary parameters and RV jitter for different datasets. Note that here we exclude the innermost planet d from the integration, assuming that it moves along a Keplerian orbit in the common system orbital plane. Its orbital period is considerably smaller than periods of other planets, as well as its mass, and therefore it does not show significant dynamical interaction with other planets on the observational timescale. Taking this planet into -body integration slows down the calculations dramatically. To ensure the gravitational influence of this planet is insignificant, we performed a similar (much longer) calculation, based on the full four-planet -body model. Almost all of the resulting best fitting parameters were practically identical, with negligible offsets of no more than of the corresponding uncertainties. The only exception was the period of this same planet d, which decreased by  day, i.e. roughly by its uncertainty. From the statistical view point, such shift should not be neglected, since it would correspond to rather large one-sigma significance level.444Throughout this paper, we will usually use the popular “-sigma” style to specify various confidence probabilities. For example, we will say that a given statement is valid at an - significance level, when corresponding confidence probability is equal to the probability for a Gaussian random variable to deviate from its mean by no more than times its standard deviation (“sigma”). The significance levels of - correspond to the confidence probabilities of, respectively, , , and (closer to means more significant). This shift is nevertheless extremely small. It does not affect the motion of the other planets at a measurable level. Since the dynamical effects due to the innermost planet are so small, we neglect them below, and integrate only the system of three remaining planets.

The apparent orbital period is also shifted due to the planetary aberration (also known as Roemer effect), which is caused by the finite light speed and is similar (in origin) to the Doppler effect (Ferraz-Mello et al., 2005). GJ876 radial velocity of  km/s (HARPS) implies that this shift should be  day (this is a “blue” shift, since the star is approaching). The cumulative correction to , due to the perturbations and Roemer effect, is  day, which is about half of the corresponding statistical uncertainty. This correction should be added postfactum to all estimations of in the fits in this paper, if such precision is necessary.

So far, we limited ourselves to the coplanar four-planet system model used by Rivera et al. (2010). When testing more complicated models, we found that the orbital fit from Table 1 shows statistically significant improvement if we add a free linear trend to the RV curve model. The estimated magnitude of this slope is about  m/(syr), which results in a quite measurable RV offset of  m/s, accumulated during the observation time span. The formal significance of this trend, as derived from the corresponding likelihood ratio statistic (Baluev, 2009), is about . Even though we have not yet took into account several important effects that we will discuss in subsequent sections, such significance is too large to be neglected without investigation. Initially, we interpreted this long-term slope as the geometrical secular RV acceleration effect, mentioned by Correia et al. (2010), which should be equal to  m/(syr) for GJ876. In fact, subtracting this predicted slope from the RV data allows to get rid of any significant trend in the fitted model. Actually, we started our research based on such corrected data, no longer bothering about any RV trends. However, A. Correia and E. Rivera later confirmed (in private communication) that both published RV time series already have this trend subtracted off. Therefore, we need to find another interpretation. We can see three equally plausible sources: another long-period unseen companion of the star, a tiny long-term instrumental drift, or some extra errors in the RV reduction pipeline. Regardless of the actual nature of this probable slope, we should try to take it into account, since it could significantly affect our results. Since its source and magnitude remains a priori unclear, we add a free linear term to the RV model. Most of our results below will refer to this model with trend, although sometimes we will make a comparison with the trend-free model. We will also return to a more rigorous estimation of the actual significance of this RV trend in further sections.

## 4 Correlated radial velocity jitter

Let us first investigate whether the currently available RV data show some residual periodicities, in addition to the current four-body model of the planetary system and possible linear trend. To do this, we utilize the common periodogram-based approach with modifications from (Baluev, 2009). The main modification is necessary to take into account the likelihood function (10), whereas the traditional periodograms (e.g. the Lomb-Scargle one) implicitly utilize the function (Baluev, 2008a). Each value of the periodogram that we are about to use, is associated with the likelihood ratio statistic, measuring how much our RV model improves, when we add a probe sinusoidal signal to it. Another important modification, which we must highlight here, is that the periodogram from (Baluev, 2009) is not the traditionally used periodogram of the fixed RV residuals to the base best fitting model. Rather, it requires a full re-fitting of the whole set of the parameters and full re-evaluation of the RV residuals for each periodogram value. This increases the sensitivity of the periodogram to faint periodicities, since we can better model the cumulative RV variation, when the probe signal indeed exists. We will refer to such periodogram as “residual periodogram”, on contrary to the “periodogram of the residuals”. The residual periodograms are obviously much more computationally-demanding than the traditional ones, especially when dealing with Newtonian -body fits.

The joint residual periodogram of all available RV data, corresponding to the base model from Table 1 shows no isolated peaks above the apparent noise level. However, it does not look like an usual white noise periodogram as well. On contrary, the noise level itself demonstrates clearly varying structure. This variability becomes especially clear, when we plot the periodogram in the linear frequency scale, instead of the logarithmic one. Fig. 2 shows such periodograms, plotted separately for the Keck and HARPS datasets (i.e., assuming the probe signal is present in only Keck or only HARPS data). The joint periodogram represents some mixture of these. We can see that in case of the Keck data there is an excessive power at low frequencies (long periods) and near the unit frequency (period close to one day), with a depression in the middle of the segment. A similar frequency distribution, although somewhat obscured by irregular variations, can be seen in the HARPS data too. These frequency spectra are very resistant with respect to various modifications in the RV curve model. The plots from Fig. 2 correspond to the circular orbit of the fourth planet, but assuming other reasonable values of and (see the next section) does not significantly affect the smoothed periodograms (although individual periodogram peaks can be affected). The influence of the long-term RV trend on these spectra looks also insignificant, as well as the influence of possible system non-coplanarity. We could not find a way to explain the non-uniform shape of the smoothed periodograms via any possible shortcomings in the RV curve model. Therefore, we need to consider another explanations.

The errors in astronomical time series are usually assumed mutually uncorrelated. Such uncorrelated sequence of errors is also known as white noise, which is called so because of its uniform frequency spectrum. A non-uniform spectrum indicates autocorrelated residuals, according to the Wiener-Khinchin theorem. This means, in particular, that the fitting methods that we used above, actually are not applicable here, since they all are based on the assumption of uncorrelated RV errors. In such a case, the correlated RV noise could be misinterpreted as some deterministic variation, and could cause unqualified systematic errors of the estimations in Table 1.

The variations of the noise level in Fig. 2 are rather large. The max/min ratios for the smoothed periodograms are and for the Keck and HARPS periodograms, respectively. Nevertheless, the apparent variations of the periodogram noise could also emerge purely by chance, because of the finite number of observations (even when the original noise is white). To demonstrate that the variations of the average noise level in Fig. 2 are statistically significant indeed, and to check in which fraction they are significant, we carry out some Monte Carlo simulations. We adopt the fit from Table 1 as the basic one, and carry out a series of bootstrap simulations of the residual periodograms plotted in Fig. 2. The bootstrap random shuffling destroys any possible correlations between different residuals, so the simulated RV noise is practically white. Therefore, the simulated periodograms should contain only natural random variations, which we should expect for the uncorrelated error noise of the same variance and distribution. Each simulated periodogram was evaluated using literally the same algorithm as real periodograms in Fig. 2, and was further smoothed to access its average level. The smoothing was done using moving average over frequency segments of  day. The result of this procedure is a bunch of simulated smoothed white noise periodograms. Then we calculate confidence ranges for the smoothed periodograms, based on the simulated periodogram set (separate confidence range for each frequency). Also, we calculate for each smoothed simulated periodogram the ratio of its maximum value over the whole available frequency range to the minimum one, and count how frequently it exceeds the same ratio of the original periodogram of the real data. It will help us to check, whether the observed variations in the basic smoothed periodograms are statistically significant or not.

The results of these simulations are plotted in Fig. 3. We can see that the simulated range of the smoothed periodograms remains almost constant in frequency. This simulated range is quite narrow for the Keck case, and more wide for the HARPS case, obviously because the number of the Keck observations considerably exceeds the number of the HARPS ones. Both periodograms of the real data do not stay in the simulated confidence ranges. In the Keck case, none of simulated smoothed periodograms could demonstrate the same or larger max/min ratio as we see in the corresponding real periodogram. This means that the noise in the Keck data is not white, with high confidence probability well above . In the HARPS case, the confidence probability is smaller – approximately  – but nevertheless is high enough to say that a similar non-whiteness probably exists in the HARPS data too, although significantly obscured by the normal random variations.

We will not try to determine here possible physical sources for such behavior of the RV noise, since this is an topic for another research. Regardless the actual sources, two main practical questions arise now: how much the noise non-whiteness affects the best fitting orbital configuration from Table 1, and how to correct this effect? These problems must be solved before we can go any further. We find that the correlated noise is already a routine issue in the exoplanetary transit searches (Pont et al., 2006). The photometric observations in these surveys often contain the so-called “red” noise component. Such name is due its power spectrum, which monotonously decreases with the frequency, thus making long-term (“red”) variations to prevail over the short-term (“blue”) ones555There is also a narrow notion of “red noise” as a synonym of the Brownian noise, that has a specific frequency spectrum . We understand the term “red noise” in a general sense, however.. Actually, the red noise is rather common phenomenon: it is not limited to only astronomy, and it is frequently faced in very different branches of the science.

The red noise indeed can easily explain the shape of the periodograms in Fig. 2: the low-frequency “hill” is just directly seen in both graphs, and another (smaller) hill near the unit frequency is its alias (caused by severe diurnal gaps in the time series). The transit observations considered in (Pont et al., 2006) allowed a simplified red noise reduction, due to their very specific distribution in time and a specific of their photometric transit models. In fact, Pont et al. (2006) could avoid, for instance, any assumptions about the shape of the noise autocorrelation function. Unfortunately, the approach used by Pont et al. (2006) is not applicable to our case. We have to model the correlated noise in full.

The details of the algorithm, which we use to eliminate the impact of the red RV jitter, are given in Appendix A. Here we only note that this algorithm is based on a certain modification of the function (10), to take the correlated noise into account. After the maximization of this new objective function, we get (at once) the estimations of the RV curve parameters (therefore, of the planetary orbital elements and masses), and several RV noise parameters: different “white” jitters for separate time series, , the common “red” jitter shared among different time series, , and the red noise correlation timescale . Such noise separation could be realistic if the red noise was actually caused by the star (e.g. by some long-living photospheric phenomena) rather than by the instruments.

## 5 Determinability of the fourth planet orbit

Further investigation of the likelihood function (10) showed that it possesses multiple maxima, which do not differ much from each other, in terms of the goodness-of-fit measures like r.m.s. or the statistic . Such phenomenon could indicate that the whole orbital model is actually ill-determined due to the lack of the RV data and/or its insufficient precision, like that was for the extrasoalar planetary system discussed in (Baluev, 2008c). However, it is not the case here. For the GJ876 planetary system, the irregularities of the likelihood function are all related to only two of the free parameters from Table 1, namely the eccentricity and the pericenter argument of the fourth planet. Fixing them at some reasonable a priori defined values makes the shape of the likelihood function very regular and practically parabolic, as it should be when the RV model is well-linearizable with respect to the remaining free parameters.

We check this for the following pairs of parameters: , , , and . The estimations of these variables possess the largest mutual correlations in the pairs, respectively , , , and (for fixed at zero). Correlated parameters should be usually more affected by non-linearity effects, as discussed in (Baluev, 2008c). For each above pair of parameters, we plot two-dimensional contours of the function (10). We choose three contours which outline confidence regions for the selected parameter pairs, corresponding to the asymptotic () confidence probabilities of . The necessary confidence probabilities were calculated from the modified likelihood ratio statistic (logarithm of the difference between the maximum and a given value of (10)), as discussed in (Baluev, 2009). After that, we compare these confidence regions with the results of numerical Monte Carlo simulations. For this goal, we use the bootstrap method as described by Marcy et al. (2005). In this method the simulated “errors” are generated by means of the random shuffling of the RV residuals to the best fitting RV curve model. This method of Monte Carlo simulation is rather widely used, because it does not require any assumptions about the shape of the error distribution, and is expected to work even for non-Gaussian errors. We introduce only two relatively cosmetic modifications to this method. First, during the shuffling we do not mess the residuals belonging to different datasets (we shuffle different datasets separately from each other). Second, after each shuffling trial, we rescale the residuals according to the total variances for the corresponding observations (which include RV jitter). We consider only the uncorrelated noise model, since random shuffling destroys any correlations anyway.

During all of the activities described in the previous paragraph, we always held and fixed at zero, to eliminate all non-linearity effects associated with these parameters. The results of these calculations are shown in Fig. 4. We can see that all of the resulting confidence domains have almost elliptic shape, and the simulated sets of points are always in very good agreement with these confidence regions. This means that when the parameters and are fixed, the remaining parameters behave almost as in the linear least-squares problem: there is a single maximum of the likelihood function (10), and this function has almost parabolic shape in the vicinity around the maximum. The estimations of the parameters possess almost Gaussian distribution, and their uncertainty estimations are pretty reliable. Notably, even the pair involving shows no significant non-linearity effects, despite rather large uncertainty in this parameter.

The behavior of the pair is severely different. The similar confidence contours for these variables look very irregular (Fig. 5) and outline two main local minima of . The first local solution is the one listed in Table 1, and the second one has larger eccentricity . This second solution appears actually a bit more likely, when we use the white noise model. But for the model with correlated noise, the first solution becomes the leading one. The picture is also sensitive to the trend term in the RV curve model. Therefore, this multi-extrema structure is model-dependent. Can we trust to this fine structure of the likelihood function at all? Most probably, we cannot: it may easily change even further, as some other non-traditional RV noise effects are discovered. The irregular details in Fig. 5 are driven by RV variations at the level of  cm/s, which is well below the internal measurement errors ( m/s). Such small variations are hardly related to the actual radial velocity of the star. They are more likely caused by unqualified systematic errors or systematic astrophysical noise in the data, which could easily have amplitude about of the random measurement noise. Therefore, the only reliable information that we can obtain here is that the eccentricity probably does not exceed , and values near are more favored than the opposite ones. All other fancy details in Fig. 5 represent just some misleading noise component, which we must try to eliminate.

This pseudo-detailed structure is not caused by some specific property of the fitting method applied here. Some points in Fig. 5 indeed yield smaller scatter of the RV residuals than the others, and this pattern is inevitably irregular. Any other classic fitting method that use the function (10) or any similar function as a single goodness-of-fit measure, would yield similar results. To understand the source of such irregular behavior of only two of the parameters, we need to trace where the information about these bad-behaving parameters comes from. To do this, we first replotted Fig. 5 for the case when the planet e is excluded from the -body RV model, and its contribution was modeled by a Keplerian function. Such plot was already pretty regular and elliptic, just as those in Fig. 4. However, the corresponding confidence domains appeared much wider than in Fig. 5, allowing values of up to . This means that the contours in Fig. 5 could not be settled by the non-sinusoidal shape of the RV variation contributed by the fourth planet, as it would be in the non-perturbed Keplerian case. Instead, the information about the likely values of and is extracted mainly from the perturbational effects, which the fourth planet imposes to the motion of other planets. Therefore, the main source of the irregularities in Fig. 5 is the dynamical interaction of the planet e with other planets in the system. This dynamical interaction allows to considerably shrink the region of admissible values of and , but by the cost of considerably irregular dependence of the RV curve on these parameters. As a result, we see some unreliable structures inside this region.

In view of the unreliable behavior of the parameters and , the best course of action for us will be to recognize that all points within some wide enough contour in Fig. 5 are equally likely. To determine, which contour should serve as a boundary, we apply a dynamical stability test to all orbital solutions spanning Fig. 5. Each point in this graph corresponds to some best fitting orbital configuration of the system. For each of these configurations, we perform the numerical integration over  yrs and check, whether a given configuration disintegrates during this test term. The integration term of  yrs is actually pretty short in comparison with cosmological timescales, but nevertheless it contains about thousand of secular periods of the system. Therefore, it roughly corresponds to  yrs, when rescaled to the subsystem of the giant planets in our Solar System, for instance. Such time segment is long enough to determine approximate stability boundaries. In Fig. 5, this stability boundary passes close to the formal confidence contour (in case of the white+red noise model). Since the stability region should somewhat shrink, when the integration time increases, we believe that it is safe to adopt the formal contour as a boundary outlining the admissible values of and .

In view of such rather poor determinability of the fourth planet eccentricity, one can ask: whether this planet exists at all? It was detected by Rivera et al. (2010) on the basis of the Keck data only, and Correia et al. (2010) did not note it in their HARPS and old Keck data. We also checked the residual periodograms for the three-planet model, and find that the fourth planet is indeed strongly supported by the Keck data, but the HARPS periodogram does not show a significant peak above the apparent noise level. One can suspect that the RV signal from this fourth planet could be actually caused by some spurious drifts in the Keck data.

However, these doubts dissolve when we look at the RV residuals to the three-planet model directly. We can see (Fig. 6), that the available HARPS data actually confirm the RV signal from the fourth planet and are in good agreement with the Keck data. This agreement just is not statistically significant yet, however it may become more significant, when more HARPS data are accumulated. At present, we have no observational basis for doubts in the existence of the planet GJ876 e.

## 6 The reference orbital fit

Using the algorithm from Appendix A, we can now obtain the full set of the parameters, taking into account the effect of the correlated noise. However, we must also address the influence of the bad-behaving parameters . Without any extra care, we will get just a few similar orbital solutions having unrealistic uncertainties. We adopt the following approach. Since other parameters behave almost linearly (see Fig. 4), we first perform a basic orbital fit, fixing the value of and at some reference realistic value (say, ). The resulting error estimation for each of the fitted parameters reflects only a fraction of the full uncertainty. We also need to estimate the uncertainty caused by the non-determinability of the parameters and . To do this, we vary and inside their admissible region (which was defined in Sect. 5) and see, how much this affects the best fitting values of other (well-behaving) parameters. This gives us the remaining part of the total parameter uncertainty (generally asymmetric).

The results of these calculations are given in Table 2 (white noise model) and Table 3 (white+red noise model). Note that when we varied and , this mainly affected the estimations of other parameters themselves, and the estimations of the statistical uncertainties remained almost constant. The first (probabilistic) uncertainties for the most of the parameters in Tables 2,3 should be quite reliable, with a very few exceptions though. These exceptions are caused, however, by the statistically unsuitable non-linear parametrization of the planetary model, rather than by the internal non-linearity of the problem itself. The estimations of bounded parameters, like the eccentricity or the RV jitter (both ) are considerably non-Gaussian and asymmetric, if the formal uncertainty range can cover the forbidden values. It is nonetheless very easy to find a better (more linear) parametrization in such cases. In particular, it is better to consider the pair instead of , as it will be done below. When the value of some RV jitter is close to zero (in comparison with its uncertainty) then it is better to deal with the squared jitter instead, with the uncertainty properly rescaled (see Table 1 by Baluev (2009)).

Looking at these orbital fits, we note, at first, that the linear RV trend appears more disputable than it seemed before: some values of and infer too low significance for (about ). Basically, the need for this trend can be eliminated by means of choosing appropriate values for and (namely, in the right and top regions in Fig. 5). On contrary, it is too early to claim that this trend does not exist at all. Other possible values of still require this RV trend. Allowing non-coplanar configurations (see Sect. 7) keeps this RV trend almost intact, so it is not easy to explain this RV trend via orbital non-coplanarity as well. We conclude that this issue can be resolved by future observations.

Second, the red noise correlation timescale is in fact poorly constrained. We can only say that has an order of days. On contrary, the magnitude of the red jitter itself, , is rather well constrained and is well separated from zero. This suggests that although the very existence of the red noise in the data looks supported, the characteristics of this jitter are still difficult to assess. Nevertheless, at present even a rough estimation of can provide important physical information.

Third, we can see some drop in the eccentricity , after we take the red noise into account. Since the apparently eccentric orbit of the planet d has an important value for, e.g., the star-planet tidal interaction theory (e.g. Ferraz-Mello et al., 2008), we should investigate the parameters and more closely. Fig. 7 contains the confidence contours for the parameters , constructed for fixed at zero, like in Fig. 4. We can see that they are perfectly elliptic and, in the white noise case, also agree with the bootstrap simulations (not shown).

Before drawing any conclusions, we need to characterize the uncertainty coming from the parameters and . Again, we adopt the contour in Fig. 5 as a safe region of admissible values for and , not paying attention to the apparent concentrations inside this region. Each point in this region refers to some orbital fit with fixed and . The values of and from these fits would mark the centers (best fit points) of the error ellipses in the plane , if we actually constructed such confidence ellipses for each admissible value. For instance, two best fit points, shown as a crosses in Fig. 7, correspond to the central points in Fig. 5. Furthermore, instead of this single point in each panel, we can construct the full set of points , which are mapped from all admissible values of . These sets are rendered in Fig. 7 as gray domains around the basic best fit points. These domains reflect how much the centers of the error ellipses in Fig. 7 may shift, while the values of and are varied inside the admissible region. Although these centers may shift, it turns out that the shape and size of the elliptic contours in Fig. 7 remain fairly constant for different and . The picture only shifts in the fairly solid way. This means that the uncertainties inferred by the bad-behaving (i.e., non-linear) parameters are well-separable from the uncertainties inferred by other (well-linearizable) parameters.

Therefore, we may treat now that the resulting uncertainty region in the plane of is constituted in a cumulative manner from the usual probabilistic confidence domains, inferred by the linear estimation theory with fixed at zero (shown as regular elliptic contours), and from the uncertainty domain inferred by the non-linear parameters and (shown as less regular small gray regions). We consider that all admissible values of and are equally possible, and therefore the gray uncertainty regions in Fig. 7 should be understood as structureless solid entities, which just indicate how the original probabilistic confidence regions should be bloated to take into account the uncertainty coming from .

We can see from Fig. 7 that is indeed inconsistent with zero, when we analyze RV data assuming traditionally that the RV noise is white. The significance of this non-zero is well above the level, even when we take into account the uncertainty inferred by the non-linear parameters and . However, it is now obvious that this apparently significant value of is likely a result of the misinterpreted red RV noise. When the red noise is taken into account, the best fitting value of moves closer to zero, and simultaneously the uncertainty regions expand. Taking into account the uncertainty coming from and , we realize that is in fact consistent with zero at the significance level of hardly above . The things remain similar when our RV curve model does not contain the linear RV trend. In that case, is non-zero at level for the white noise model, but for the correlated model this significance drops to the same level. Although we still cannot retract the values of as large as , we nevertheless have no observational evidences that is actually non-zero. Further observations can eventually solve this question for sure, but at present it is too early to claim that the non-zero eccentricity was confirmed. The apparent non-zero value of reported in previous works represents basically an effect of misinterpreted red noise in the RV data.

## 7 System non-coplanarity

One of the primary goals of this paper was to characterize the mutual non-coplanarity between planets b and c or at least to put some limit on it. In the non-coplanar case, we have two separate variables for the osculating inclinations and , and also a pair of extra parameters determining the orientation of their ascending nodes, and . Two latter parameters, however, cannot be estimated independently, because the problem is invariable with respect to arbitrary rotation around the line of sight. In fact, we can determine only the difference . The main quantity that we are interested in, in view of the orbit non-coplanarity, is the mutual orbital inclination for the planets b and c, which is a function of , , and .

We must note that during this work we will deal, most probably, with small values of , comparable to its statistical uncertainty. Such situation reveals many practical difficulties, because our estimations of will have significantly non-Gaussian distribution. Similar troubles arise for planetary eccentricities, when they are small enough (comparable to their uncertainties). In the both situations, the issue is rather formal, however. The non-Gaussian behavior of such estimations is caused by the bad choice of the parametrization, rather than is inferred by the RV model itself. In case of the eccentricity, the problem is usually eliminated if we consider, instead of the eccentricity, the pair of parameters . Estimations of these parameters usually have almost Gaussian bivariate distribution, even when is small. All apparent problems in this case are caused by the trivial singularity of the polar coordinate system in the point . It is very likely that the mutual inclination should be affected by a similar singularity at , which could be eliminated by means of transition to the variables like , where is an extra auxiliary variable, determining some extra orientation angle related to the mutual orbital inclination. Note that the variables and should be independent in the sense that for any pair of and there should exist a valid orbital configuration. Due to that requirement, we cannot choose to be equal to, e.g., the angle between the planetary ascending nodes . For any , the value of cannot exceed , and such non-constant constraint will imply nothing except extra difficulties. We could use some orientation angle in the Laplace plane of the system. Such choice would be physically justified and also symmetric with respect to the two planets involved, but it would mess the geometric parameters (like inclinations) with planetary masses, which is not desirable.

We adopt the following purely geometric definition of the auxiliary angle . Given two abstract “first” and “second” planets, we can determine two mutual orbital nodes of these planets. We choose the node in which the second planet ascends over the orbital plane of the first planet, and define as the orientation angle of this node in the plane of the first orbit, counting it from the usual ascending node of the first orbit. This definition is schematically illustrated in Fig. 8. With a help of the classical spherical trigonometry, it is not hard to derive the formulae expressing the new osculating angles via the original parameters . They look like:

 sinIcosΦ = −sini1cosi2+cosi1sini2cosΔΩ, sinIsinΦ = sini2sinΔΩ, cosI = cosi1cosi2+sini1sini2cosΔΩ. (12)

The inverse transition is given by the equalities:

 sini2cosΔΩ = sini1cosI+cosi1sinIcosΦ, sini2sinΔΩ = sinIsinΦ, cosi2 = cosi1cosI−sini1sinIcosΦ. (13)

Obviously, it is impossible to express all three original orientation parameters via only and . We must also know to obtain and . Therefore, such definition of is not symmetric with respect to the two planets: the first orbit serves as a reference plane. In case of GJ876, we choose the planet b to be this “first” planet, since its RV amplitude is considerably larger, and thus its orbital plane orientation is determined with better precision. The planet c will be the “second” planet, and the orbit of the planet e is assumed to lie in the common Laplace plane of the system (at the epoch of osculation). The planet d was also assumed to move in the Laplace plane, though this assumption only affects the mass estimation of this planet.666To make these definitions more rigorous, we must also mention that in the Jacobi coordinate system, that we adopt here, different osculating orbits have different reference points, so they are no longer confocal, and this is not reflected in Fig. 8. These small displacements do not have practical significance, however.

Given the formulae (12,13), we can easily carry out constrained RV fits with and fixed at any desired values. Technically, such constrained fitting may be done, for instance, by means of expressing and via , , and (thus eliminating and from the set of free parameters). Therefore, we can perform a series of such constrained fits on some regular grid of and and then to plot the resulting likelihood contours in the plane, e.g., . Such plot basically visualizes the confidence regions for these variables, similar to those regions that we have already constructed in Fig. 4.

Results of these calculations are shown in Fig. 9 for both noise models (white and white+red). It is notable that in the white noise case the inclination shows some non-zero value at the significance level of , although this significance becomes somewhat smaller when the uncertainty in is included. These non-coplanarity signs are further softened, when the red noise is taken into account. In this case, is consistent with zero at level, even when the and uncertainties are neglected. Eventually, we conclude that available RV data for GJ876 are fully consistent with the coplanar configuration (). Nevertheless, we can obtain some informative upper limit on the possible non-coplanarity from Fig. 9: the angle likely cannot exceed . However, because of the significant prolateness of the error ellipses in Fig. 9, the upper limit on significantly depends on the value of , i.e. on the orientation of the orbital nodes for planets b and c. So large values of as can be accepted only if the corresponding orbital planes intersect each other rather close to the sky tangent plane (i.e., close to or ). When the intersection nodes are far from the sky plane, the mutual inclination is unlikely to exceed .

From the non-coplanar three-planet fit by Correia et al. (2010) we find and from their Table 2 and , from their Table 3. Rivera et al. (2010) only mentioned that their non-coplanar four-planet fit yields . Both groups agree that there is no significant mutual inclination between the planets b and c. Our results do not contradict to such estimations. We have to admit, however, that both works discussed the non-coplanarity issue very briefly, and they omit exact uncertainties for . Correia et al. (2010) gave some uncertainties for (about ), but did not supply the necessary correlations, disabling us to derive the inferred uncertainties for and/or .

Previously, Bean & Seifahrt (2009) also tried to determine the orbit non-coplanarity between the planets b and c, based on the old Keck RV data from (Rivera et al., 2005) and HST astrometry data from (Benedict et al., 2002). Their non-coplanar orbital fit corresponds to and in our notation. These values also agree with confidence regions in Fig. 9. We must note that although Bean & Seifahrt (2009) utilized the astrometic measurements, they, however, took into account neither the annual systematic errors in the old Keck data (Fig. 1), nor the very existence of the fourth planet, nor the red noise in RV data. They also do not mention whether they removed the secular acceleration effect from the RV data they used. Therefore, their results cannot be directly compared with ours. Due to this, the possible effects from the astrometry data by Benedict et al. (2002) remain not fully clear. We do not expect, however, that these effects are large. Indeed, note that the estimations of, for instance, the absolute node longitudes from (Bean & Seifahrt, 2009) possess rather large uncertainties of . Since the astrometry data are responsible for these large uncertainties almost exclusively, we may suppose that such data should not constrain very much the values of , and probably most of the non-coplanarity information now comes from radial velocities anyway. The role of the astrometric data was nevertheless more important in (Bean & Seifahrt, 2009), since they used old RV data, which allowed much larger of . We do not use the astrometric data here, because their error properties are in fact poorly assessed. Bean & Seifahrt (2009) noted that the actual scattering of the astrometric residuals is considerably different from the stated instrumental uncertainties. Also, these data can easily contain a correlated component. We think these effects are difficult to estimate reliably, due to the relatively small size of the astrometric dataset.

In the first, purely white, case reflected in Fig. 9, the results of the bootstrap simulations (not shown) are in good agreement with the confidence regions plotted, just like in Fig. 4. In the second, white+red, case, the bootstrap simulations are not helpful, since they destroy any noise correlation effects anyway. However, it is very likely that the formal statistical reliablity of the confidence regions in this case should be so high as for the white noise model. The only remaining question is how well the particular noise model, used in the algorithm from the Appendix A, allows to eliminate the effects coming from the RV noise correlateness. To check this, we replotted the confidence contours from Fig. 9 assuming a bit with different noise models. First, we checked the picture remains practically the same for different simple noise correlation functions (, , and ). After that, we probed different splitting of the red part of the RV noise between the Keck and HARPS data. We checked the cases when the HARPS noise is purely white, and when it has its own red component, not tied to the Keck one. In these cases, the changes in the confidence regions are larger, roughly similar to the difference between left and right panels in Fig. 9. This may indicate that the effect of the red RV noise still may be not taken into account in full. It seems from the current data, that the HARPS red jitter may be a bit smaller than the Keck one. We cannot use a more accurate noise model, however. The Keck red RV jitter alone looks pretty estimatable without the HARPS data, but the HARPS red jitter is, on contrary, poorly separable. We have to live with that problem until more HARPS data are acquired. Anyway, all noise models tested so far, do not imply a significant non-coplanarity of the system and place almost the same limits on it.

## 8 Long-term dynamics

Now we can investigate the long-term dynamical evolution of the planetary system. We choose two best fitting orbital configurations, corresponding to the osculating (solution I) and (solution II). Both configurations are obtained assuming the white+red RV noise model. We tracked the evolution of each orbital configuration over  Myr term (the innermost planet d was not taken into account during the integration). The relative energy error was about in the first integration, and about in the second one. The energy error was walking around these values all the time and did not show a notable accumulation effect. This is exactly what we should expect from the symplectic integrator. We must note, however, that our dynamical analysis here is still rather preliminary and in future it is better to perform the integration over longer terms and to take into account the short-period planet d using, e.g., an averaged Hamiltonian method (Farago et al., 2009).

Both orbital configurations appeared stable during the integration time. The evolution of the parameters related to the massive planets b and c is fairly regular and is close to the apsidal corotation resonance state. The main secular period of the system – the period of the pericenters revolution – is close to  yrs in the both cases. The perturbations from the fourth planet are rather small, although they add some minor chaotic component in the motion of the massive planets. The long-term evolution of this fourth planet itself looks, on contrary, considerably chaotic, in terms of its orbital eccentricity at least. This eccentricity, however, remained bounded from the upper side by (solution I) and (solution II).

We also considered the evolution of the following resonant variables, corresponding to the individual two-planet resonances:

 scb,c=2ψb−ψc−ωc,scb,b=2ψb−ψc−ωb, sbe,b=2ψe−ψb−ωb,sbe,e=2ψe−ψb−ωe, sce,c=(4ψe−ψc)/3−ωc,sce,e=(4ψe−ψc)/3−ωe. (14)

These definitions of resonant angles are derived from the general definitions from (Beaugé et al., 2003). Since we consider only coplanar configurations here, the angles and are replaceable by and (as we explained in Sect. 2.1). The first pair in (14) corresponds to the 2/1 MMR between the planets b and c, the second pair – to the same resonance between b and e, and the third pair – to the 4/1 MMR between c and e. All eleven two-planet resonant angles studied by Rivera et al. (2010) can be expressed via the six variables (14). Namely,

 φcb,c=−scb,c,φcb,b=−scb,b,φcb=scb,b−scb,c, φbe,b=−sbe,b,φbe,e=−sbe,e,φbe=sbe,b−sbe,e, φce0=−3sce,c,