Enhanced models for stellar Doppler noise reveal hints of a 13year activity cycle of 55 Cancri
Abstract
We consider the impact of Doppler noise models on the statistical robustness of the exoplanetary radialvelocity fits. We show that the traditional model of the Doppler noise with an additive jitter can generate large nonlinearity effects, decreasing the reliability of the fit, especially in the cases when a correleated Doppler noise is involved. We introduce a regularization of the additive noise model that can gracefully eliminate its singularities together with the associated nonlinearity effects.
We apply this approach to Doppler timeseries data of several exoplanetary systems. It demonstrates that our new regularized noise model yields orbital fits that have either increased or at least the same statistical robustness, in comparison with the simple additive jitter. Various statistical uncertainties in the parametric estimations are often reduced, while planet detection significance is often increased.
Concerning the 55 Cnc fiveplanet system, we show that its Doppler data contain significant correlated (“red”) noise. Its correlation timescale is in the range from days to months, and its magnitude is much larger than the effect of the planetary body perturbations in the radial velocity (these perturbations thus appear undetectable). Characteristics of the red noise depend on the spectrograph/observatory, and also show a cyclic time variation in phase with the public Ca II H&K and photometry measurements. We interpret this modulation as a hint of the longterm activity cycle of 55 Cnc, similar to the Solar 11year cycle. We estimate the 55 Cnc activity period by yrs, with the nearest minimum presumably expected in 2014 or 2015.
keywords:
techniques: radial velocities  planetary systems  methods: data analysis  methods: statistical  stars: activity  stars: individual: 55 Cnc1 Introduction
After the Doppler discovery of the famous planet 51 Pegasi b in 1995 (Mayor & Queloz, 1995), a large progress was made in the exoplanets detection field, and the number of robustly detected exoplanet candidates have recently exceeded (Schneider, 1995). However, the development in this field have encountered a natural barrier that is difficult to overcome. This barrier is related to the stellar activity: even if we decrease the instrumental errors of the radialvelocity (hereafter RV) measurements to a level of cm/s, we would not usually obtain the r.m.s. below m/s due to the intrinsic astrophysical instabiltiy of even quiet Solartype stars. Only the quitest stars demonstrate the Doppler stability at a level below m/s.
Any random noise can be suppressed, to a certain degree, by accumulating a large number of observations. However, this path is increasingly hard, because the data analysis get affected by an increasing number of various subtle noise effects that must be therefore taken into account. In modern works on the analysis of exoplanetary Doppler data we can see that the authors are in a fierce fight with each such systematic effect, hoping to reduce the statistical uncertainties and detection thresholds, and to eradicate various noise misinterpretation errors (e.g. Dumusque et al., 2012; Baluev, 2013c; Hatzes, 2013; Nelson et al., 2014; Tuomi et al., 2014).
One of such effects is the correlated Doppler noise discovered recently in the precision RV data for some stars, having a correlation timescale of d in average (Baluev, 2011; Tuomi & Jenkins, 2012; Baluev, 2013c; Feroz & Hobson, 2014; Tuomi et al., 2014). This correlated, or ‘red’, RV noise have not got a good theoretic discussion yet, and depending on the case it may be inspired by various instrumental effects as well as by stellar activity phenomena (e.g. star spots).
In addition to the detection of previously unknown noise components, other issues related to the noise models appear. These issues are related to the constructions of robust models for the known noise components. Unfortanutely, so far we have no good physical understanding of the structure of the Doppler noise at the accuracy level of m/s. Models that are currently used to approximate the white or red RV noise often represent only rather primitive formulae that were chosen mainly on the basis of mathematical simplicity and convinience with the use of only minimal or no knowledge of the underlying physics. Although some sources of the RV jitter like the stellar astroseismology oscillations were identified and investigated long ago (e.g. O’Toole et al., 2008), the topic of the RV jitter structure still represents a rather speculative and suggestive matter with insuffiecient experimental/observational base, so the difficulties arising with accurate modelling of the RV jitter are objectively defined.
Acknowledging this, we consider the issue of mathematical usefulness of the RV noise models. Definitely there are multiple models that can approximate the required general physical behaviour equally well, but have considerably different mathematical properties. Which one should be selected in practice? In the case of exoplanetary Doppler data we should take a special care of the statistical robustness of the derived RV fits. Under the fit robustness we understand here the degree of nonlinearity of the fitting task. Clearly, if we have two models that both adequately reproduce the desired physical behaviour, we should select the one that generates a smaller nonlinearity effects, because unjustified nonlinearity in the fitted model usually only triggers various troubles.
What undesired consequences a nonlinear model can impose? First of all, it is a violation of the wellknown asymptotic (largesample) properties of the classic maximumlikelihood estimation methods: the parametric estimations become considerably nongaussian (sometimes even multimodal with comparable modes), the likelihoodratio test and all statistical tools that rely on it become not applicable, and the maximumlikelihood estimations become biased and loose their optimality (the asymptotic property of maximum efficiency). Also, highly nonlinear models require more iterations of the fitting algorithm, in comparison with the welllinearizable cases. This results in an increased demand of computational resources.
One can object that modern statistical dataanalysis methods like the Bayesian inference (Ford, 2006; Gregory, 2007a, b; Tuomi & AngladaEscudé, 2013; Nelson et al., 2014) can gracefully handle even highly nonlinear models, contrary to the classic maximumlikelihood approach. While recognizing that the Bayesian methods indeed are currently the bestdeveloped statistical tools for solving nonlinear estimation tasks, we emphasize the practical difference between correct accounting for nonlinearity, which is what Bayesian methods can do, and suppressing the nonlinearity on the other hand. The ability of adequate handling of nonlinear models does not mean that the nonlinearity effects are removed. Thus, even with the Bayesian approach we should try to use, whenever possible, a welllinearizable model as a preferred alternative to a highly nonlinear one.
The construction of a welllinearizable and simultaneously physically adequate model of Doppler noise is thus the primary goal of the paper. In Sect. 2 we explain how a noise model can generate any nonlinearity of the fit, how the traditional model with an additive RV jitter may become highly nonlinear, and what is a welllinearizable noise model (only considering the white noise case). In Sect. 3 we describe the new “regularized” model of the white Doppler noise that eliminates the weaknesses of the additive model. In Sect. 4 we consider the regularity issues raised by the red noise models. In Sect. 5 we discuss several issues concerning various periodograms, and how they are affected by nonlinearity. We also give there slightly extended results concerning the periodogram significance levels for models involving the red noise. In Sect. 6 we present testcase Monte Carlo simulations demonstrating that the regularized noise model may lead to a dramatic improvement of the fit robustness in some practical cases. In Sect. 7 we present the results of the new analysis of the Doppler data for the star 55 Cancri that hosts a fiveplanet system (using the new noise modelling technique). In Sect. 8 we discuss the hints of a decade activity cycle revealed in these Doppler data and the relationship of these results with the observations from other longterm monitoring projects that involved 55 Cnc.
2 White noise models and nonlinearity
Let us assume that we deal with the time series containing measurements (e.g. radial velocities), acquired at the times , and having the random errors . The model that should be fitted with is designated as . We now consider several models of the noise and the consequences they can carry for the fit of .
In this section we assume that represent white Gaussian noise, i.e. that they are independent normally distributed variates. In this case to set a ‘model’ of the noise means to express the variances through some noise parameters by means of a relationship
(1) 
where the function represents the noise model that we are speaking of. This framework was already considered in general terms, e.g. by (Baluev, 2009), where it was suggested to fit the curve parameters and the noise parameters simultaneously by means of the maximumlikelihood approach. This can be achieved through the maximization of the loglikelihood function:
(2) 
Here we consider a few special cases of the noise model , depending on a single scalar parameter .
The classic models of the noise assume that either have known standard deviations or at least known statistical weights . In the first case we have
(3)  
(4)  
(5) 
Now the task is reduced to the classic leastsquare fitting as .
In the second case we deal with the arbitrarily scalable noise:
(6)  
(7) 
where is an unknown parameter and is given by (5). We also call this classic noise model as ‘multiplicative’ one. Clearly, now the likelihood function can be maximized over and separately, so that we have , and is again obtained through the leastsquares fitting.^{1}^{1}1The reader may notice that this estimation of is actually biased, while the classic unbiased expression should be , where is the number of the curve parameters, . The necessary bias correction is performed in (Baluev, 2009) by modifying the likelihood function (2), but we do not highligt this change here, keeping more simple formulae for demonstration.
In this paper we are mainly concerned by the effects of the nonlinearity of the task. A strictly linear curve model , where represent some known functional base, would imply that is expressed though a quadratic form, so that its minimization becomes trivial. If the model is formally nonlinear, but the function still can be approximated by a positivedefinite quadratic form, we say that the leastsquare task is welllinearizable. If cannot be approximated by a quadratic function, especially if it possesses multiple local minima, we say that the associated leastsquare task is essentially nonlinear. However, it is not trivial to classify the noise models in the same manner. Here we should base on the likelihood function (2) rather than on the noise model itself. Namely, we should say that the noise model is ‘linear’ if the loglikelihood is expressed by a quadratic form in terms of the noise parameter . The model is not supposed to be linear here in the strict mathematical sense.
Now let us return to the classic multiplicative model (6). Is it linear in this sense or not? Apparently, even if we consider fixed, the relevant loglikelihood (7) is not quadratic in . However, we may note that this loglikelihood always has only a single maximum in , and it tends to both for and . Therefore, for a fixed we can easily find a parametric replace , such that the loglikelihood function is quadratic in . We cannot yet say that the inferred likelihood function is strictly quadratic in both and , because the replace also depends on . However, in practice we may consider the values of in a small () vicinity around the best fitting solution, so we may say that (7) can be rather accurately, although not exactly, represented by some quadratic function in . Thus, the multiplicative noise model (6) can generate little nonlinearity effects. The cumulative nonlinearity of the maximumlikelihood task is thus dominated by the nonlinearity coming from the curve model .
After that, let us consider the following ‘additive’ noise model from (Baluev, 2009):
(8)  
(9) 
where designate some known uncertainties, while the unknown parameter is responsible for the additional ‘jitter’. The model (8) was originally constructed assuming that physically represents the variance of some additional noise component, implying (Wright, 2005). However, here we do not require to be necessarily positive. The model (8) may remain quite meaningful even for . The negative value of means that the real errors are systematically smaller than the stated uncertainties , because e.g. these were overestimated due to some imperfections of the spectrum analysis pipeline. In practice such cases exist indeed (Baluev, 2009), and therefore the ability of the noise model (8) to handle such cases is pretty desirable, even though its original motivation becomes nonphysical for .
Nonetheless, not every negative value of is meaningful. We may notice an undesired irregularity in the likelihood function near , which is the minimum admissible value of . When approaches to this value, one of tends to zero, but contrary to the multiplicative model, the others remain separated from zero. This implies that the relative weight of the single observation with the minumum increases infinitely, so this observation plays an increasingly dominant role. This is obviously an undesired behaviour: we have no any reasonable justification to infinitely increase the contrast between the weights. This singularity does not always show itself up. Normally, the residual corresponding to the observation with minimum does not vanish, and the values of near this limit are naturally penalized by the term in (2). This term tends to at a greater speed than the other term, , tends to . The resulting loglikelihood thus tends to , indicating an unsuitable solution. However, we obtain a degeneracy when this residul vanish, leaving only the logarithmic term in the game. In this case the entire loglikelihood function would tend to , indicating an absolutely (!) preferrable solution.
This singularity effect is illustrated in Fig. 1. As we can see, the singularity is very sensitive with respect to the residual , so to be actually trapped, the likelihood maximization algorithm should meet in its path a small enough residual simultaneously with a negative value of . In fact this singularity effect was already mentioned in (Baluev, 2009), where also a workaround solution to the issue was suggested. The workaround was to cut off the singularities by truncating the additive model (8), replacing it with
(10) 
and using the general likelihood (2). This forbids the maximumlikelihood fitter to travel to the domain where any of becomes smaller than . For example, the truncation with is applied internally in PlanetPack (Baluev, 2013a), i.e. it is not allowed to reduce the uncertainties more than by the factor of . In practice this pretty mild condition allows the fitter to be never trapped in the singularity, although it is still possible to have the best fitting value of lying at the very boundary of the truncation.
Although the truncated model (10) solves the numerical issues of the fitting by cutting off the likelihood singularites themselves, it does not eliminate the severe nonlinearity effects that expand over a wider domain around the singularity. The observations with smaller values of still receive a disproportionally increasing relative weights for without good physical reasoning. The other observations cannot have significant effect on the fit in such a case. This basically means that the effective number of observations involved in the game gets decreased, pumping up the nonlinearity of the task. However, we do not want to just ban the negative values completely, because as we disscussed above, negative allows to handle practical cases with overestimated errors. Besides, any abrupt cutoffs of the parameters introduce other nonlinearity effects, so the model (10) is still just an ugly workaround rather than a neat solution to the issue.
To suppress the nonlinearity effects generated by the additive noise model for we should construnct a new ‘regularized’ noise model with better mathematical properties in this domain. In the domain this new model should behave like the additive model (8), while for it should behave similarly to the multiplicative model (6), keeping the weight ratios of the observations constant.
The graph of the additive model (8) consists of two straight rays: for and for . Given this, we define our regularized model as a hyperbola constructed on these rays as on asymptotas. Mathematically, such regularized model may be expressed as
(11) 
with the likelihood given by (2). Here is some dimensional scale factor, is a fittable adimensional jitter parameter, and is a small regularization parameter. For the model (11) is strictly equivalent to the additive jitter model. For nonzero and the model (11) is asymptotically equivalent to the additive model, with the relative error decreasing as . For we obtain , which becomes equivalent to the multiplicative model. Also, for (11) we have , meaning that our regularized model for coincides exactly with the additive model for (regardless of the values of ).
3 Regularized model of the additive jitter
The regularization (11) is not a unique choice, although this is one of the simplest mathematical functions that suits our needs. Note that the original multiplicative (6) and additive (8) models can be represented in the same form as in (11), substituting the following functions in place of :
(12) 
We can see that in comparison with , the function offers a smooth nonbreaking decrease of the variance for , which removes any pecularities from the likelihood function and thus may reduce the effects of nonlinearity. Note that the relative contrast between the weights of the observations keeps bounded by the contrast between . Thus we also eliminated the pathological property of the additive model that forced only one or a few of the RV measurements to dominantly contribute in the final fit if .
The regularized model (11) has two finecontrol parameters, and . They are not independed: the replacement preserves the model. Therefore, we may set to any reasonable dimensional value and then select the nondimensional parameter . We set equal to the geometric mean of , and then after some experimenting we find that good compromise values for are located in the range from to . The graphs of the regularized model (11) for a few values of are plotted in Fig. 2.
The jitter parameter in (11) and (12) is just a dummy adimensional variable. For the multiplicative and additive jitter models, we define the relevant dimensional RV jitter using the formula
(13) 
where is the parameter of (12). Here we follow the convention that negative corresponds to negative (implying overestimated ). The same definition could be used for the regularized model (11) too. However, in this case we may want to make the regularized and the additive jitter models interchangeable in the sense that switching between them should not change the implied jitter estimation too much. With (13) this requirement is not satisfied. For example, as we noted above, the regularized model for coincides with the additive one for . In this case we would expect to have , but we obtain from (13) a nonzero value , which is not always satisfactory. For smaller negative , the bias in is further increased, e.g. the case should correspond to the effective jitter of about , rather than to minus infinity implied by (13).
Let us consider this effect in more details. From (11) it follows that
(14) 
Now let us assume that actually obey the additive model (12) for some , and we try to approximate these using a slightly different model (11). That is, we seek a replacement such that . From (14) follows that for each individually such a representation can be made exact if we substitute in (11)
(15) 
However, we should find a single common value of approximating all simultaneously. Such value of can be expressed using the same formula (15), but replacing all by some their mean value . Moreover, if we choose wisely, this mean value should be close to unit. Therefore, replacing in (15) by unit, we obtain a suitable ‘best fitting’ representation for :
(16) 
Thus the entire procedure of determining the jitter estimation looks as follows. First of all, we perform the maximumlikelihood fit of the regularized model (11) via the parameter . Then we derive from (16) the matching value of :
(17) 
which represents the argument of an additive model that approximates the fitted regularized one. Finally, we substitute this value in place of into the definition (13). When is defined using this procedure, its global minimum value is , which maps to . The value maps to , and this now correctly corresponds to . For large positive this corrective offset becomes negligible.
This also implies that the actuall modelling errors of (11) are smaller than it may seem from Fig. 2, because in this figure no correction was applied to the argument . Considering the regularized model in (11) as a function of , like , we would obtain exact coincidence with the additive noise model for . For or , and with , we obtain that the modelling errors of remain below for (see Fig. 3). Therefore, even is still a small enough value that allows (11) to efficiently mimic the additive model in the domain . However, for the divergence grows quickly, allowing the regularized model to have much better mathematical properties. These properties lead to an improved robustness of various statistical estimations and to a better predictability of various statistical thresholds.
4 Regularity of the red noise models
The activityrelated Doppler variations in solarlike stars were investigated relatively long ago (e.g. Saar & Donahue, 1997; Saar et al., 1998), but the available RV data amount and quality disabled any detailed modelling of this Doppler noise until recent years. In the most of the exoplanetary RVanalysis works in the past, the Doppler “jitter” was modelled as a white noise. Now it is clear that a well detectable correlated (or “red”) noise is often present in the exoplanetary Doppler data. The number of known stars with significant rednoise component in the RV data is growing (Baluev, 2011, 2013c; Tuomi & Jenkins, 2012; Feroz & Hobson, 2014; Tuomi et al., 2014), and these examples demonstrate the typical noise correlation timescales of d in average. Apparently, this red noise is excitated by the stellar activity effects (Dumusque et al., 2012; Hatzes, 2013; Robertson et al., 2014).
In (Baluev, 2011, 2013c) and (Feroz & Hobson, 2014) the red noise in Doppler data was taken into account by means of adopting a realistic (although rather arbitrarily selected) parametrized model for its covariance structure. The following twocomponent model of the correlated RV noise was applied:
(18) 
where , , and are free fittable parameters, while is the adopted shape of the correlation function. We used , corresponding to the socalled “red” noise.
The likelihood of the correlated data is different from (2). Approximating the red noise by a stationary Gaussian random process with the covariance matrix (18), the likelihood can be expressed by
(19) 
with the residuals vector same as in (2). This likelihood is a function of three noise parameters , , and , and of the usual parameters of the RV curve, . The joint estimation of all these parameters can be obtained by finding the maximum of (19), as usually.
We can see that the noise model (18) is potentially vulnerable with respect to the irregularities discussed above, because we did not put any constraints on or , except for an implicit condition that the resulting noise covariance matrix V should be positive definite. Note that with the model (18) the cases of become very frequent, because a significant fraction of the available residulas scatter may be consumed by the “red” term with . Moreover, because we now have two noise components, we may face two types of singularities. They appear in the domains or . A singularity in the likelihood (19) formally occures when passes zero. The singularity of the white component can be removed by the truncated model (10) or preferably by the new regularized model (11).
What concerns the red noise component, here the solution is even more simple: we should just forbid to attain negative values, as they are nonphysical. In this case we do not have any justification to possibly allow , like we had with . The constraint can be easily implemented e.g. by replacing and treating as a primary fittable parameter.
5 Likelihoodratio periodograms
Periodograms represent an extremely important analysis tool for Doppler exoplanets detections. Besides, periodograms are rather sensitive to the nonlinearity of the adopted models, see e.g. (Baluev & Beaugé, 2014). In fact, we originally suspected the nonlinearity issues of the additive model (8) because of the too frequent disagreement between the simulated periodogram significance levels and their expected analytic approximations.
The general likelihoodratio periodogram was formally defined in (Baluev, 2009), and it basically represents the difference between maximum values of loglikelihood obtained for two RV curve models, with and without a probe sinusoidal signal. This is a general definition valid both for various noise models, including white or correlated.
We would like to highlight an important internal property of these periodograms that often makes them superior over those periodograms that are still used in the majority of other exoplanetary dataanalysis works. This property becomes important when we deal with multiplanetary fits, e.g. when we have already detected one planet, and need to verify whether the data contain yet another signal. In this case we should compare the base (or null) model, involving only the signal of the known planet, and the alternative model involving the signals of the known planet + additional sinusoid. The widespread solution to this issue is to fit the base model, compute the residuals of this best fit, and then treat these residuals as a new standalone time series, applying to it the same periodogram that was used to detect the first planet (i.e., assuming a zero base model). It was clearly demonstrated by (AngladaEscudé & Tuomi, 2012) that this approach is inoptimal and leads to a decreased detection power. The second periodogram here does not take into account the possibility that the first best fit may be inadequate or inaccurate due to the contribution of the putative second planet. Thus, the residuals obtained after the first fit may be corrupted by the inaccuracies of the base model. A better solution is to adaptively adjust the parameters of the first planet in the process of computing the second periodogram. In other words, the second periodgram should involve the fit of the original data (instead of the singleplanet residuals) with the entire twoplanet model (i.e. first planet + putative sinusoid, rather than just the sinusoid), comparing this fit with the fit of only the base model. This modification leads us to the likelihoodratio statistic and the associated likelihoodratio periodograms from (Baluev, 2009).
We call such periodograms as “residual periodograms”, contrary to the usual “periodograms of residual”, while AngladaEscudé & Tuomi (2012) calls them “recursive periodograms”. The value of the residual periodogram is always greater than or equal to the value of the periodogram of the fixed residuals, because the first one involves a complete maximization of the likelihood function, while the second one corresponds to another, nonmaximum, value of the likelihood. Simultaneously, the noise level of the residual periodograms remains approximately the same, making the residual periodograms more sensitive to faint signals. This was practically demonstrated by AngladaEscudé & Tuomi (2012).
Note that when the base model is expressed by just a constant the same issue leads to the construction of the “floatingmean periodogram” (FerrazMello, 1981; Cumming et al., 1999; Zechmeister & Kürster, 2009). The floatingmean periodogram is a more efficient tool than the Lomb (1976)Scargle (1982) periodogram of a centred time series. Basically, the residual likelihoodratio periodograms just represents a further generalization of the floatingmean periodogram to the cases with more complicated models of the RV curve and nontrivial models of the RV noise.
The likelihoodratio periodograms also allow an extension of another type, useful when dealing with heterogeneous data, e.g. RV data for the same star coming from independent teams. In this case we may be interested in detecting possible descrepancies between the different data sets. This can be achieved by constructing periodograms referring – in some sense – to individual data sets. We construct such datasetspecific periodograms by comparing the base model and the alternative model, obtained by adding a sinusoidal signal to the base model only for data points that belong to a particular subset. For other data points the alternative model is set equal to the base one. The resulting likelihood ratio statistic allows to reveal variations belonging to only a particular dataset, still using the full statistical power of the entire time series to fit the base model. This method was introduced in (Baluev, 2011) and its practical efficiency was further demonstrated in (Baluev & Beaugé, 2014). Clearly, this approach is more efficient then e.g. a plain restriction of the analysis to only a single RV data set.
The crucial matter in the periodogram analysis is the determination of the statistical significance of the periodogram peaks that we reveal. The analytic theory for the periodogram significance levels is given in (Baluev, 2008). Summarizing, the false alarm probability associated with a given maximum peak, can be approximated as
(20) 
where is the periodogram maximum and is the adopted frequency range to scan. Also, is an effective timespan of the data, which is proportional to the weighted variance of :
(21) 
When dealing with periodograms that are tied to a particular data set in a heterogeneous time series (as explained above), we should just restrict the summations in (21) to the relevant subset.
The sign ’’ in (20) means that the function in (20) represents simultaneously an upper bound and an approximation for the . The latter approximation is asymptotic: for large (small ) its error decreases. It is important that the formula (20) was originally constructed only for strictly linear models with only the frequency of the probe sinusoidal signal allowed to be nonlinear. However, for nonlinear but linearizable models the approximation (20) often should work in the asymptotic sense, assuming (Baluev, 2009). Therefore, possible deviations between the simulated and analytic curves, especially the ones that break the inequality in (20) may indicate the presence of a significant nonlinearity. See also additional discussion in (Baluev, 2014).
Likelihoodratio periodograms can be easily defined for models involving red noise. The significance levels for such periodograms can be constructed using the same approach. Similarly to the white noise case, we may construct the quadratic Taylor decomposition of the loglikelihood function (19) (near its maximum). Such a decomposition allows to transform all the formulae to a shape similar to the classic linear regression task considered in (Baluev, 2008). In the end, it appears that the formula (20) is again asymptotically applicable even for models with correlated (but still well linearizable) noise, with a slightly generalized value of . Namely, we should now substitute in (21)
(22) 
where V is the compound noise covariance matrix of the RV data (corresponding to the best fit), is a vector containing units, and is a vector containing all . In the rednoise case we should also satisfy the extra requirement with being the noise correlation timescale. Finally, when dealing with periodograms referring to individual RV subsets of the heterogeneous time series (see above), we should restrict the vectors and in (22) to this particular subset, but the matrix V should remain intact.
The values of usually appear somewhat larger for the red noise, but it seems that the difference in the predicted levels often appears negligible (see e.g. Fig. 7 discussed below).
6 Testcase simulations
We performed some Monte Carlo simulations revealing the effect of the jitter model on various statistical distributions associated to the best fitting models. In this section we consider relatively demonstrative examples involving RV data for several wellknown exoplanetary systems.
Originally we noticed the effects of the nonlinearity in the additive jitter model when considered the simulated distributions of periodograms. For large , such periodograms should asymptotically obey the LombScargle periodogram distribution derived in (Baluev, 2008), provided that the RV curve model is welllinearizable. But contrary to our expectations, the simulations show rather large deviations too often, even for very large data sets and rather simple RV models. As we demonstrate below, in many cases such deviations are owed to the nonlinearity of the additive jitter model, rather than to the ‘genuine’ nonlinearity of the RV model. We show that our regularized model of the jitter can significantly improve the statistical behaviour of such periodograms.
Additionally to periodograms, we noticed that the additive noise model may significantly exacerbate the nonlinearity impact in the cases that involve a correlated Doppler noise. We demonstrate below that our new regularized noise model may be very helpful in such cases too.
6.1 Doppler data for 51 Pegasi
In the first testcase we consider the public radial velocity measurements of 51 Pegasi: ELODIE ones (Naef et al., 2004) and Lick ones (Butler et al., 2006). We first set a base model of these data. This includes a Keplerian almost sinusoidal variation induced by the planet 51 Peg b (Mayor & Queloz, 1995), the longterm linear trend (Butler et al., 2006; Baluev, 2009), and a sinusoidal annual variation in the ELODIE data responsible for some instrumental drifts or data reduction errors (Baluev, 2009). For the sake of symmetry between Lick and ELODIE, we added a sinusoidal anuual variation to the model of the Lick data too, although contrary to ELODIE its amplitude was not statistically significant. In these annual variations the frequency was considered as a fixed parameter, while the trigonometric coefficients were freely fittable. The compound base model looks like
(23) 
where year (fixed), and the last term of the model (on a separate line) is the usual Keplerian signal due to 51 Peg b. The index in (23) is equal to either or , depending on whether the model is computed for an ELODIE or Lick data point. Basically this model contains only a couple of sinusoids (the Keplerian signal has almost zero eccentricity and thus is almost sinusoidal) plus a linear trend. Meanwhile, the number of the measurements in each of the datasets is very large, so the RV curve model should be pretty linearizable. We would not expect that so simple model can generate any nonlinear overfit effects with so large time series. With this RV curve model, we consider three different models of the Doppler noise: the classic multiplicative one (6), the legacy additive ones (8,10), and the new regularized one (11).
Apparently no pitfalls were encountered with the fitting of any of these noise models. The best fitting values of the RV curve parameters are pretty expected and unremarkable, so we do not detail them here. We only note that while the derived ELODIE best fitting jitter value of m/s was well separated from zero, the Lick RV jitter appeared consistent with zero within its uncertainty: (m/s). The singularity of the additive noise model is located at the jitter value of . In the Lick data this minimum uncertainty is m/s, corresponding to only a sigma separation between the derived jitter estimation and the singularity. The average uncertainty value is m/s. Therefore, we may expect that the Lick data represent, due to their negative jitter, a hidden source of the nonlinearity and nonrobustness of the fit, even if the RV curve model itself looks well linearizable. On contrary, the ELODIE data should not produce any such problems. For them we have m/s with almost all being equal to or m/s (minimum is m/s).
Note that here we applied the modification to the likelihood function (2) that allows to perform a preventive bias reduction in the derived jitter values (Baluev, 2009). Therefore, it is not the statistical bias that makes the estimation of the Lick jitter close to zero and even negative. Actually, thanks to the large number of the data, this bias correction appears rather small, so the results would be qualitatively the same even if we used the classic likelihood function (2).
Now we consider various likelihoodratio residual periodograms of these RV data. First of all, let us look at the Lick periodogram (Fig. 4). In these plots we show two selected chunks in the period axis, containing the false peaks caused by the irregularity of the additive model. We can see that the plain additive model (8) fails here, and the fitting algorithm is trapped in the singularity. The heights of the bogus peaks should be actually infinite, but the fitting was stopped due to some of the internal stopping criteria (likely, maximum number of iterations reached). This periodogram contains many such bogus peaks, spanning a wide period range, and basically rendering the periodogram itself useless. This behaviour is largerly corrected by the truncated additive model (10) that effectively cuts the false peaks.
However, the periodogram for truncated additive model is still not always smooth in the vicinities of the problematic periods. It is clear that although the singularities themselves were cut off, the nonlinearity that always resides in the domain around a singularity is still not eliminated and not suppressed. Although this nonlinearity does not generate any obvious undesired effects in this individual periodogram, we may suspect that it may introduce certain biases when dealing with large randomized periodogram ensembles. This means that various probabilistic estimations, like statistical uncertainties and significance thresholds may get disturbed. For example, the average noise level on the periodogram may be increased in comparison with what we would expect from a welllinearizable case, resulting in a decreased detection efficiency.
Finally, our regularized noise model (11) allows to obtain an entirely smooth periodogram. Because we deal here with negative jitter (i.e. overestimated instrumental errors), this periodogram stays close to the periodogram that was obtained for the multiplicative model (6). On contrary, if the best fitting jitter was positive, the periodogram with regularized noise would be closer to the one with additive noise. Therefore, our regularized noise model is able to automatically decide whether the additive model is realistic or not for a particular data, and smoothly select the suitable behaviour.
Now let us proceed to periodogram significance levels. For each of the three noise models mentioned above we perform a Monte Carlo simulation of the Lick, ELODIE, and joint periodograms in the frequency range from to day. Based on each of these simulations we construct the empirical cumulative distribution function (CDF) of the periodogram maximum, and compare this CDF with its analytic approximation (20). The results are shown in Fig. 5. In these graphs we plot the ratio as a function of , for mentioned periodograms. Since Monte Carlo simulations infer some random uncertainties that should be honoured, we also plot (in grayscale) the sigma confidence ranges for this ratio, computed by means of the weighted KolmogorovSmirnov test (see Appendix A).
First of all, we can see that the ELODIE periodograms (first panel) do not show any signs of nonlinearity with any of the noise models. In fact, these plots only reveal that the formula (20) behaves exactly as it was supposed to behave in the linear case: it represents an asymptotic approximation to the when , simultaneously serving as an upper limit for larger s. Also, we may note that all three curves generated by different noise models almost coincide with each other within the Monte Carlo uncertainties. This is because of the small spread of for the ELODIE RV data: in fact, almost all of these are equal to either m/s or m/s. In this case we have all in (11,12) close to unit, implying that all three noise models are practically equivalent to each other. In this case the weights of the observations become almost equal for any of these models, and we just need to assess the unknown common noise variance.
In the second panel of Fig. 5 the results for the Lick periodograms are shown. In this case the curve for the truncated additive model does not obey the analytic model (20). It is remarkable and rather disappointing that even the inequality in (20) got broken: the simulated may be even twice as large as . This indicates that Lick periodogram with additive noise possesses an increased noise level, in comparison with what we might expect from a linearizable fit. As we suspected, the truncated modification (10) is unable to suppress this nonlinearity. However, our regularized model is clearly able to suppress it, making the analytic formula (20) useable again. The results for the periodogram of the joint ELODIE+Lick data (third frame in Fig. 5) look similar to the Lick case.
The summary is that the additive RV noise model, including its truncated version, often generates an excessive noise on the periodograms in the cases when the estimated jitter value is negative or consistent with zero (taking into account it uncertainty). Obviously, this noise excess may only decrease the detection efficiency, and increase the statistical uncertainty of the usual RV curve parameters. This noise excess is generated by the nonlinearity of the RV noise model only, and thus it can appear even when the RV curve model is very simple and does not show any significant nonlinearity in itself.
Fortunately, this nonlinearity can be efficiently suppressed by replacing the additive noise model with the regularized one. Moreover, with the regularized noise model we noticed up to per cent speedup of the computations. The likely explanation is that to fit a welllinearizable model a smaller number of LevenbergMarquardt iterations is needed than to fit a highly nonlinear model (in the ultimate limit of strict linearity even a single iteration would be enough).
6.2 Doppler data for GJ 581
The RV data for this star contain a significant correlated noise component with a correlation timescale of d (Baluev, 2013c; Tuomi & Jenkins, 2012). Here we revist the results of (Baluev, 2013c) in view of the improved noise models introduced above.
In (Baluev, 2013c) we used the approach described in Sect. 4 above, but it was unintentionally allowed to have negative, unless it breaks the positive definiteness of the noise covariance matrix V. Also, we used the plain (not even truncated) additive model of the white noise in this case. As we now demonstrate, these issues caused significant corrupting effect on some of the Monte Carlo simulations from (Baluev, 2013c). Although the Monte Carlo trials generating were seldom, they triggered drastically large values of the maximized likelihood function, introducing significant distortions in various simulated distributions.
This is shown in Fig. 6, where we compare simulated distributions of the likelihoodratio statistic for different noise models. The curves labelled with “1” are for the plain model (18) and were taken from (Baluev, 2013c). The curves labelled with “2” correspond to the model where the white noise component is expressed by the truncated model (10) with , and the red noise is also expressed by a “truncated” model, in which negative are disallowed. The curves labelled with “3” and “4” differ from the second case only in the model of the white noise component, which is set to the regularized with or to the multiplicative one, respectively. In the left panel of Fig. 6 we show the simulations done for the twoplanet RV curve model, and in the right panel we include three planets (this is the maximum number of planets that can be detected in this system with only the Keck RV data). We can see that in both the cases, the curve “1” is very different from the others, indicating that this simulation was significantly affected by the nonlinearity that emerges for negative or .
Similarly to the 51 Pegasi testcase, this nonlinearity results in reduced significance estimations and increased statistical uncertainties. Fortunately, it can be suppressed by choosing another noise model. Therefore, this is in fact a “fake” nonlinearity caused by the bad mathematical behaviour of an incarefully selected noise model. Here we do not spot any remarkable difference between the truncated additive, regularized, and multiplicative models of the white noise component, although the difference with the plain model (18) is obvious. A detailed investigation revealed that in the twoplanet case both types of the singularity with and with may appear in the simulations (although not simultaneously of course), while for the threeplanet model the Monte Carlo trials with do not occur. We may conclude that in practice any of the two noise terms in (18) can raise irregularity issues.
Since we now have an efficient method to suppress the nonlinearity of these rednoise fits, we used it to recalculate the significance of the fourth planet in the system, GJ 581 d. In (Baluev, 2013c) this significance was estimated by a rather marginal value of sigma at most. Given the simulation results discussed above, we may expect that this significance may be underestimated due to the poor quality of the RV noise model, and in true this planet can be confirmed with a greater confidence. We rerun the periodogram simulations from (Baluev, 2013c) assuming the regularized model (11) for the white component and constraining the red jitter parameter to only positive values. The correlation shape function was adopted exponential again.
origin of the estimation  only HARPS data  joint HARPS+Keck data,  joint HARPS+Keck data, 

“shared” red noise  “separated” red noise  
Monte Carlo estimations from (Baluev, 2013c)  sigma  sigma  sigma 
Monte Carlo estimations from the present work  sigma  sigma  roughly sigma 
Analytic formulae (20,21,22)  sigma  sigma  sigma 
Normal quantiles corresponding to the relevant periodogram s are given. The “shared” and “separated” red noise models correspond to the cases when the parameters for the Keck and HARPS red noise are bound with each other or are fitted independently, see (Baluev, 2013c) for details. These results do not take into account the recent Robertson et al. (2014) work.
We summarize our results in Table 1. We can see that now the detection significance levels of the planet d are remarkably increased. Moreover, the new Monte Carlo estimations are now in a much better agreement with the analytic approximations that assume a linearizable task. The simulated significances are slightly larger than the analytic ones, formally breaking the inequality (20), still indicating the presence of some residual nonlinearity in the models. However, the deviation is now modest and pretty negligible. This indicates that the detection of the GJ 581 d signal is now statistically robust.
Also, our updated results do not call back the conclusion of Baluev (2013c) that adding the correlated noise component drastically reduces the periodogram peak corresponding to the planet d. In fact, after applying our regularized model the height of this peak in the periodograms of the real data changes negligibly in comparison with what we had in (Baluev, 2013c). All the distortions are living in some fraction of the simulated periodograms, increasing the periodogram noise. Also, the conclusions of Baluev (2013c) concerning the more speculative planet candidates f and g remain intact: the relevant periodogram peaks just disappear after turning on the RV noise correlations in the model. Here we further confirm that these planets likely do not exist or at least not supported by the current RV data.
Although our new reanalysis of the GJ 581 d signal suggests that it appears formally more significant than in (Baluev, 2013c), we still do not insist on its existence. During the refereeing of the present paper, a new work by Robertson et al. (2014) appeared, where the authors claimed that this signal owes to the stellar activity manifesting itself in the Doppler data for GJ 581. However, they demonstrated that the actual correlated noise in the GJ 581 HARPS data is not stationary, likely showing some cyclic variability in time. Contrary to this, our analysis, both in (Baluev, 2013c) and here, assumes a constant activity level and constant red noise characteristics. If the actual red noise is varying, its reduction in the constantlevel approximation might appear incomplete. This means that the significance estimates in Table 1 should decrease back when the variations of the red noise characteristics are taken into account. It seems that GJ 581 is in this regard reminiscent of the 55 Cnc case considered below. However, in this paper we had no aim to consider the GJ 581 data in so deep details. Here we only used its data as a testcase to demonstrate the nonlinearity of the white and red noise models and verify that our regularization scheme indeed allows to suppress most of this nonlinearity.
6.3 Doppler data for HD 82943
Of course, our regularized noise model is not a magic stick that would suppress the nonlinearity everywhere it is applied. We must remember that in some cases the nonlinearity of the fit genuinely belongs to the model of the RV curve, rather than of the RV noise. Such an example is offered by the HD 82943 system (Tan et al., 2013; Baluev & Beaugé, 2014). The regularized noise model was designed to handle more gracefully the cases in which the estimated jitter is consistent with zero or even negative. However, in the HD 82943 all jitter estimations (see Baluev & Beaugé, 2014) are positive and well separated from zero. For example, Baluev & Beaugé (2014) present the simulation of the periodogram of the Keck data for this star (see their Fig. 12). This simulation shows significant deviation from the analytic prediction (20), owing to the fit nonlinearity. Now we recomputed this simulation based on the regularized noise model (11), but no visible change was revealed. In this case, the nonlinearity is caused by the traditional sources: the complexity of the RV curve model with a relatively small number of the Keck RV observations.
7 Correlated Doppler noise in the 55 Cancri RV data
Nelson et al. (2014) presented a thorough investigation of the updated RV data for this star. This investigation included the latest Keck (Nelson et al., 2014), Lick (Fischer et al., 2014), HET (HobbyEberly Telescope, Endl et al. 2012), and HJST (Harlan J. Smith Telescope, Endl et al. (2012)) measurements. Although many subtle effects were considered in the analysis by (Nelson et al., 2014), e.g. a correlation of Doppler noise over a single night due to stellar oscillations, they did not address the possibility of RV noise correlations at the timescale of a few days or more. They only stated (without an explanation) that the assumption of uncorrelated RV errors “is usually an excellent approximation when observations are separated by one or more days”. But as we discussed in Sect. 4, this assumption may appear wrong, because at the precision level of m/s we frequently deal with RV correlations on the average time scale of d.
Here we undertake an attempt to address this issue of Doppler noise correlations over the timescales greater than d, based on our new approach involving regularized noise models. We perform a reanalysis of the abovementioned public RV measurements from Keck, Lick, HET, and HJST, also adding to them the ELODIE dataset from (Naef et al., 2004), which Nelson et al. (2014) did not include. Following to Nelson et al. (2014); Fischer et al. (2014), we separate the Lick dataset into several subsets, acquired using different Dewar, according to the information given in the Fischer et al. (2014) Lick RV database. We need to note that the Lick dataset used by Nelson et al. (2014) is a bit larger than the one published Fischer et al. (2014). It seems that the difference is due to the treatment of the RV measurement series that were obtained in a single night. Apparently, Fischer et al. (2014) replace such series by some average value, because these data do not contain so close , while Nelson et al. (2014) deal the full (unpublished) data, treating the correlations inside the series close data points in two distinct ways (their “Case 1” and “Case 2”). Here we only have access to the public Fischer et al. (2014) data. From the discussion by Nelson et al. (2014) we learned some additional information about the Fischer et al. (2014) RV data. The starting group of RV measurements that are labelled with Dewar code 8 is actually related to an “early Doppler era” and should be treated separately from the ending group with formally the same Dewar code. Also, it appears that these earlyera data should be separated in two even smaller subsets ( points) that refer to different Dewar and may have different statistical characteristics. However, due to the small number of these observations we cannot model their jitter separately: such an attempt frequently leads to parametric degeneracies. We therefore decided to share a single jitter model between these two RV subsets, still allowing them to have different RV offsets. We must note that these RV measurements anyway have a small effect of the entire RV fit due to their small number and relatively large uncertainties.
At first, we considered the unperturbed (multiKeplerian) model of the RV curve, involving all known planets orbiting 55 Cnc, and approximating the RV noise in the abovediscussed datasets by the regularized model (11). These periodograms based on the regularized white noise model are shown in the left column of Fig. 7. For each of these periodogram we also plot a triplet of , , and sigma noise levels computed using the analytic formula (20).
We can see that they look very differently for different datasets. In the most cases we can see a clearly nonuniform frequency distribution, indicating that the Doppler noise in not white.^{2}^{2}2With the logarithmic scale of the period axis used in Fig. 7 the white noise should appear monotonically growing to the shortperiod end of the axis. If such a pattern is not fulfilled, the noise is not white. Only the Lick Dewar 6 subset and the Keck1 subset are consistent with the white noise. Simultaneously they agree well with the predicted noise levels. The other RV subsets show remarkable excess power at the periods d and in a narrow period range around d. This is a typical picture that the socalled red noise generates: the power excess at long periods is aliased to the vicinity of the d period due to diurnal gaps existing in the RV data. Simultaneously, these periodograms do not comply with the predicted noise levels expected for the white noise. We can see many peaks that rise well above even the threesigma level, but these peaks are different for different datasets, indicating some noise effects rather than an actual RV periodicity.
We modelled this red noise separately for each RV dataset, excluding the Lick6 and Keck1 ones (they preserved only white components). The model of the red part of the noise was analogous to the one used for the GJ 581 case (see above), constraining . The white part noise was modelled by the regularized model again. The crosscorrelations between different RV subsets were set to zero. The periodograms based on such compound noise model are plotted in the right column of Fig. 7. We can see that the red noise is largerly suppressed in these periodograms. In the most cases the periodograms become consistent with the predicted noise levels, indicating an almost complete reduction of the red noise. Only the Keck2 and ELODIE data still contain some remaining nonwhiteness and marginally significant peaks.
To verify that the nonlinearity of the rednoise model is small enough and our statistical methods a legal here, we perform Monte Carlo simulations of the likelihoodratio statistic. These are done similarly to the GJ 581 case above, and the results are shown in Fig. 8. We can see that in the whitenoise case the agreement between the simulated and asymptotic distribution is very good. For the rednoise model only some minor deviation appears. Therefore, we believe that both the RV curve and RV noise models do not generate significant nonlinearity in the fit, and our results should be pretty reliable in terms of their statistical robustness.
We also consider the impact of planetplanet Newtonian perturbations on the fit. We consider best fits based on either Keplerian or Newtonian (coplanar edgeon) model of the RV curve, and calculate the difference between the inferred RV curve models. Additionally, we compute the difference between the best fitting RV curves for the white and red noise models (neglecting Newtonian perturbations). These differences are plotted in Fig. 9. We can see that Newtonian perturbations have typical value of only m/s over the observation time span, significantly smaller than the RV uncertainties. Moreover, we can see that the impact of the red noise exceeds the Newtonian perturbations roughly by the factor of . This means that the RV data available for 55 Cnc are still unable to reveal any signs of Newtonian perturbations in this planetary system. Therefore, these Newtonian perturbations have only a negligible effect on the periodograms and simulations discussed above. Moreover, it may even appear that we wrongly interprent some phantom rednoise variation as a hint of Newtonian perturbations. From the other side, (Nelson et al., 2014) show that the fit accuracy is enough to have a significant difference between the osculating and averaged orbital elements, exceeding the estimation uncertainty by the factor of a few. These conclusions do not contradict to each other: the presence of detectable offsets between the osculating and averaged orbital elements does not necessarily imply the presence of detectable perturbations in the RV curve (except for a secular linear drift in the planetary mean longitudes that act like a period biasing).
planetary orbital parameters and masses  

planet b  planet c  planet d  planet e  planet f  
[d]  
[m/s]  
[]  
[]  
[]  (fixed)  
[]  
[AU]  
parameters of the data sets  
Lick early 1  Lick early 2  Lick 13  Lick 6  Lick 8  
[d]  
[d]  
[m/s]  
white [m/s]  
red [m/s]        
[d]        
[m/s]  
r.m.s. [m/s]  
Keck 1  Keck 2  ELODIE  HET  HJST  
[d]  
[d]  
[m/s]  
white [m/s]  
red [m/s]    
[d]    
[m/s]  
r.m.s. [m/s]  
general characteristics of the data and fit  
[d]  
[d]  
[m/s]  
The parameters have the following meaning: orbital period , RV semiamplitude , eccentricity , pericenter argument , mean longitude , derived planet mass and semimajor axis . The osculating orbital elements are in the Jacobi reference frame described in (Baluev, 2011), but the innermost planet was involved neither in the Nbody integration, nor in the construction of the Jacobi reference frame. The values of and were derived assuming the common orbital inclination of and the mass of the star . The uncertainty of was not included in the uncertainties of the derived values. The parameter is the constant RV offset, and is the estimated magnitude of the white or red RV jitter (the “Lick pre1” and “Lick pre2” datasets share the same jitter). The values represent the correlation timescale of the red noise. The white jitter is approximated by the regularized model with the specified jitter scale factor (see text). The quantities and represent a dataset usual timespan and the effective timespan needed to compute periodogram noise levels (see text). The goodness of the fit is tied to the modified likelihood function as explained in (Baluev, 2009). The integers and are the number of observations and of the free parameters, respectively.
However, considering all pro and contras, we prefer to deal with the osculating elements and Newtonian model, and we give this final 55 Cnc fit in Table 2. This fit involves the correlated RV noise model and Newtonian coplanar edgeon model of the RV curve. Comparing this fit with the fit by Nelson et al. (2014), we do not spot any significant difference, except for the planet f eccentricity, . We obtain , while Nelson et al. (2014) give an almost zero value. This change appears due to the rednoise model. More detailed investigation shows that this is still consistent with zero: we obtain only sigma deviation according to the relevant likelihoodratio statistic. In fact, this eccentricity appears very poorly determined, offering only an upper limit of .
The orbital inclination in Table 2 is fixed at (edgeon) because otherwise it appears undeterminable: we obtain that the value of the maximized likelihood function changes rather little down to as small as (see Fig. 10). However, (Nelson et al., 2014) give the estimation with the uncertainty of . We view this uncertainty too optimistic, and we guess that it may owe e.g. to the prior distribution adopted for by Nelson et al. (2014) in their Bayesian simulation. Note that the planet 55 Cnc e is transiting its star, so its orbital inclination to the sky plane, , is necessarily close to . In Fig. 10 we vary the common orbital inclination of the remaining planets , thus allowing the entire system to be noncoplanar. According to the results by Nelson et al. (2014), the fiveplanet system becomes unstable when the mutual inclination between the orbital plane of 55 Cnc e and the common orbital plane of the other planets is in the range from to . As far as is close to , the inclination , considered in Fig. 10, cannot be smaller than .
Finally, we note that in (Baluev, 2013b) we have done a brief period analysis on the 55 Cnc Lick RV data from (Fischer et al., 2008), claiming an uncertain detection of an additional d signal. The period of d was not confirmed by Nelson et al. (2014) with updated data, and now we do not confirm it too. This d signal now looks like an artifact of the red noise that was ignored in (Baluev, 2013b).
8 Evidences of an activity cycle in Doppler data of 55 Cancri
It is remarkable that according to our fit, the ELODIE and HJST datasets possess considerably different noise correlation timescale d, whereas the other datasets have d. This may indicate that the correlated noise in the ELODIE and HJST subsets has essentially different physical nature. Perhaps, it might have an instrumental origin in ELODIE and HJST, while in the other data it may be caused by the star activity.
In fact, there is no unique explanation of the differences in the correlation parameters between RV datasets. First of all, we do not really know how much of the red noise was generated by the instrumental sources and by the star. But even if the red noise is entirely caused by the star activity, its numerical parameters may largerly depend on the characteristics of the spectrograph or telescope, or on the spectrum analysis software. Also, the differences in the correlation parameters may indicate some temporal variability of the star activity level.
In the case of 55 Cnc the data are very heterogeneous, coming from different observatories and teams, so it is difficult to say which effect is responsible for the difference in the rednoise characteristics. However, if we leave only the Lick and Keck data that were obtained by means of more or less similar observation technique and spectrum analysis methods, it appears that the red noise level was close to a minimum in Lick6 and Keck1 observations, while in Lick13, Lick8 and Keck2 data it was near a maximum. This infers an activity cycle of yrs, with a single observed minimum near JD2452500 (about 2002.5) and two observed maxima near JD2450000 (about 1996.0) and JD2455000 (about 2009.5).
Various monitoring observations available in public literature on 55 Cnc reveal clear hints of the longterm activity variation in this star. Results by Henry et al. (2000) suggest a cyclic variation in the Ca II H+K activity indicator with a single observed maximum in 1995. They do not cover the complete cycle with good confidence, but the activity minimum can be suspected some time after the end of their observations timespan in 1999. A longterm downtrend was noted by Marcy et al. (2002) in the Strömgren photometry between the dates JD2450500 and JD2452500. A later work by Fischer et al. (2008) presents an updated plot of year photometry series. In this plot we do not see any longterm flux variation, because it was subtracted by the authors. However, some longterm cyclic variation in the magnitude of the photometic scatter can be easily noted. This scatter is likely tied to the spotting activity of the star, which is therefore variable. If this interpretation is correct, we can see some excessive activity starting from the beginning of the data and continuing until JD2451400. It is followed by a depressed activity epoch lasting until JD2453200. After that, an increased activity regime returns, continuing till the end of the dataset (JD2454200). The apparent minimum is likely located near JD2452000, and the second maximum is likely near JD2454000. Unfortunately, Fischer et al. (2008) do not give any period estimation or even period analysis of the longterm flux variability, although they note that “55 Cnc clearly exhibits yeartoyear variations in mean brightness with an amplitude of 0.001 mag over a timescale of several years or more”.
These results clearly provide an independent support to our hypothesis that the Doppler red noise in the 55 Cnc Lick and Keck RV data owes to the star’s magnetic activity cycle, similar to the Solar year activity variation. Likely, this red noise is caused by spots evolving on the stellar surface, and the correlation timescale of d may reflect the intrinsic stability of these spots or spot groups. This timescale is also in a rough agreement with the typical lifetime that we have for sunspots or sunspot groups (days to weeks). Also, this timescale can be related to star rotation. The rotation period of 55 Cnc is d (Fischer et al., 2008), and after d the star would make of a single revolution, causing significant changes to the spot pattern seen on its visible disk, and consequently to the measured radial velocity.
Given this considerations, we undertook an attempt to estimate the period and phase of the longterm variability of the red noise component in Doppler data. To do this we replace the noise model (18) with the following modulated modification:
(24) 
Here and are the period and the phase of the modulation, and now corresponds to the maximum rednoise magnitude. Note that when we have the variance of the red noise component proportional the modulation factor of . This assumes that in the epoch of minimum activity the red noise vanishes entirely. It is easy to verify that the matrix in (24) keeps its positive definiteness if it was positive definite in (18). In (24) we assume a regularized model of the white noise, , and require that .
We include in the analysis only the Lick and Keck data, except for the very old and rather inaccurate “earlyera” Lick data points. The ELODIE and HJST data have too different value of that may indicate a different source of red noise correlations, possibly having instrumental nature. Besides, the ELODIE data are rather old and frequently reveal a systematic annual variations possibly owing to some reduction errors (Baluev, 2009; Baluev & Beaugé, 2014). Presumably, HET data are located near the epoch of minimum activity, but in their periodogram we can see a clear power excess at longer periods (Fig. 7). However, it is mainly generated by a single high peak at the period of d. The time span of these data is very small (the same d), and this variation may be due to some systematic instrumental or datareduction drift. Previous studies revealed a systematic annual variation in the HET data for HD 74156 (Baluev, 2009; Wittenmyer et al., 2009; Meschiari et al., 2011). Whether the same is true for the HET data from (Endl et al., 2012) or not, this RV subset remains too suspicious, and we belive it is unreliable for this very subtle work.
Thus, we are left with highquality RV subsets: Lick13, Lick6, Lick8, Keck1, Keck2. In the model (24) we assume that the whitenoise term may have different value of for different RV subsets. The values of and are assumed the same across the entire RV data, while and are assumed the same for all Lick data and for all Keck data, but are allowed do differ between Lick and Keck (this looks reasonable based on the data from Table 2). The cross correlations between different datasets is still set to zero, because otherwise we have to deal with a nontrivial question, what and we should set for these crosscorrelation matrix blocks of V, in particular bearing in mind the requirement of positive definiteness of V. If the red noise is indeed caused by the star activity then this assumption of zero cross correlations is probably not very physical, but in practice different Lick and Keck RV subsets either do not overlap or overlap little, so this assumption does not make significant difference to the final results. The RV curve is approximated by a fiveplanet Keplerian model. Bearing in mind the results presented in Fig. 9, the effect of the Newtonian planetplanet interaction is much smaller than the effect of the red noise, and can be therefore neglected.
Lick 13  Lick 6  Lick 8  Keck 1  Keck 2  
white [m/s]  
red [m/s]  
[d]  
r.m.s. [m/s]  
[d]  
[]  
general characteristics of the fit  
[m/s]  
At first we performed the maximumlikelihood fit of the described RV model. We obtained the best fitting values of and with no remarkable changes in the planetary parameters. We give the other RV noise estimations in Table 3, omitting the planetary system parameters. However, it appears that this period estimation is uncertain. The Doppler data cover the minimum activity epoch of 2001 well, but the two maximum activity epochs of 1995 and 2009 are covered only partially. Therefore, with only the Doppler data the uncertainty of becomes very asymmetric without a reliable upper limit. This is demonstrated in Fig. 12, where we show the likelihood regions for the parameters . The likelohood contour levels shown there formally correspond to the asymptotic 1,2,3sigma significance levels of the relevant likelihoodratio test.
For comparison, we also plot over this graph a set of Monte Carlo points. They were obtained by assuming that the best fitting model of the real RV data with the implied noise covarince matrix are true, then simulating in accrodance with this covariance matrix a set of synthetic RV datasets, and fitting each of these datasets with the same model. If the model was significantly nonlinear, the simulated points would deviate from the contours of the likelihood function. However, instead of this we can see a good agreement. This likely indicated that the model is well linearizable, the fit is statistically robust, and from the statistical point of view we can safely rely on the uncertainties derived with the likelihoodratio test. Unfortunately this does not yet mean that the best fitting estimations and uncertainties of and are indeed so reliable in the practical sense. They appear severely modeldependent and until we understand well the underlying mechanism of how the red noise is generated in the Doppler data, we should treat these estimations and confidence domains only as a rough suggestive reference.
In Fig. 12 we show similar likelihood regions for two other equivalent pairs of parameters: (minimum activity epoch + 1st maximum activity epoch) and (minimum activity epoch + 2nd maximum activity epoch). The activity minima and maxima are understood here as the minima and maxima of the Doppler red noise sinusoidal modulation. Anyway, we can see that only the lower period limit is more or less useful, while the upper one is very vague: even the periods as large as yrs are allowed.
However, the Ca II H&K data from (Henry et al., 2000) yield the epoch of the first activity maximum with rather good accuracy, and this can be used to dramatically improve our constraints on . In the left frame of Fig. 12 we show the horizontal hashed band that refers to the year of 1995, which is the maximum activity year from (Henry et al., 2000). We refitted our model with a constraint that the phase of the modulation is zero at JD2449800 (some date in 1995). This suggested a refined estimation yrs (uncertainties are from the likelihoodratio test). Thus, the 55 Cnc activity cycle is similar to the Solar cycle or slightly longer. The next year of 2015 is a predicted year of minimum activity of 55 Cnc, although we emphasize that this prediction, as well as the period estimation above, are sensitive to model assumptions, in particular to the assumption that the activity is strictly periodic.
In Table 3 we can see that the Doppler red noise in the Lick and Keck data have considerably different parameters: both the correlation timescale and the maximum rednoise variance are significantly larger for the Lick data. The difference in looks rather natural, because the Lick/Hamilton spectrograph likely can “see” only the bigger, and thus more longliving, spots and spot groups on the star disk due to the weaker techinical characteristics. However, the difference in the red noise magnitude is rather paradoxical. It is unexpected that a weaker instrument is more sensitive to a subtle spot activity. We have several explanations to this:

Lick data contain an additional rednoise component that may have instrumental nature.

The star was more active in the maximum of 1995 than in that of 2009.

Some hidden technicalities can suppress the red noise in the Keck data. This may be broadly analogous to the longexposure averaging effect on the shortperiod astroseismology harmonics.
The estimated activity period is suspiciously close to the period of the outermost and most massive planet d. Can this modulated red RV noise be caused by e.g. Newtonian perturbations from this planet on the motion of the inner planets? According to the results of Fig. 9, the effect of the correlated noise exceeds the Newtonian perturbations for by a factor of . Thus these two effects may become comparable only for as small as , scaling the real planetary masses up enough. However, even in this case there are several arguments against such interpretation of the red RV noise in the data:

The dynamics of the system on so short time span is definitely regular, implying that the frequency spectrum of the Newtonian perturbations in the radial velocity should contain only a few isolated frequencies tied to the planetary orbital mean motions and their linear combinations, or to the secular periods. The red noise that we can see in periodograms of Fig. 7 does not demonstrate this.

The magnitude of perturbations and other their characteristics should be the same for the Lick and Keck data, but this does not holds for the red noise that we see.

The inclination of is inconsistent with the results by McArthur et al. (2004), who give for the outermost planet based on HST astrometry + RV data available that time. The magnitude of the astrometric displacement should scale roughly as (due to an increase of the real planetary masses), meaning that if was only then McArthur et al. (2004) would be wrong by a factor of , or even larger due to the projection effect. This does not seem likely, even if the results by McArthur et al. (2004) were based on an old and inaccurate value of .

There is a clear correlation of the derived Doppler activity cycle with the Ca II H&K and photometry observations discussed above.
Therefore, we treat the coincidence between the activity cycle and the period of the outermost planet as accidental. Interestingly, the Solar activity cycle is close to the Jupiter orbital period too, and in the 55 Cnc system the period of the planet b is close to the star rotation period. From the other side, some discussion is available in the literature investigating the relationship between the Solar activity and planetary orbital motion (Abreu et al., 2012; Cameron & Schüssler, 2013; Tan & Cheng, 2013). Our estimation of the 55 Cnc activity cycle may offer some fresh source material to a work of such type.
9 Conclusions
The radialvelocity method of exoplanets detection has encountered a natural limit set by the intrinsic activity of even relatively quiet Solartype stars. Recent studies reveal that careful modelling of the Doppler noise may be a key to the future progress of the radialvelocity exoplanet detections. For example, Dumusque et al. (2012) were able to reveal the planet orbiting Cen B only after suppressing the star’s rotation and longterm activity effects in their Doppler data (although this detection currently appears modeldependent, see Hatzes (2013)). In (Baluev, 2013c) the planet GJ 581 e could be confirmed in the Keck data only after taking the red noise into account. From the other side, an application of advanced noise models may lead to planet retractions too: see the cases of GJ 581 f and g in (Baluev, 2013c; Tuomi & Jenkins, 2012), and of GJ 581 d in (Robertson et al., 2014), and likely the case of the GJ 667 C multiplanet system in (Feroz & Hobson, 2014).
The present paper further emphasizes the importance of a wise selection of the Doppler noise model and the need to increase the complexity of the latter. For example, in Table 2, about 2/3 of the model parameters is related to the Doppler noise, while only 1/3 of them refers to the exoplanetary system. Moreove, this work reveals that it is not enough to think only of the formal accuracy and physical justification of the RV noise model. Additionally to this, we must take care of mathematical smoothness of the resulting likelihood function, and carefully define the admissible parametric domain, paying a special attention to the regions where the noise model is formally not physical but still mathematically defined.
The 55 Cnc case appears even more plentiful in what concerns the usefulness of nonstandard noise models. Other solartype stars exist, in which an activity cycle was detected by means of the Doppler RV monitoring (Dumusque et al., 2011), but 55 Cancri is likely among the most detailedly investigated ones. Probably, more practical methods stellar activity monitoring are the the Ca II H&K lines measurements (Henry et al., 2000; Robertson et al., 2013) and highaccuracy photometry (Fischer et al., 2008). Since these methods allow more easy and accurate way to track the stellar activity, they may be useful in construction of a more reliable model of the red Doppler noise. Instead of performing a pure Doppler fit that would rely on some particual temporal model of the longterm modulation, we may try to subtract this modulated red noise based on its correlation with less expensive observations. This would be generally similar to what Dumusque et al. (2012); Tuomi et al. (2014) suggest, but the activity index is correlated with the magnitude of the red RV noise instead of the raw radial velocity. Thus it is still necessary to answer the question, which approach is more appropriate in practice. In any case, a good reduction of the stellar activity hints in Doppler data would likely increase the accuracy and reliability of the derived exoplanetary model for the stars like 55 Cnc.
An RV fitting algorithm based on our regularized jitter model is now implemented in the PlanetPack package (Baluev, 2013a) and will be available in the forthcoming version PlanetPack 2.0. This version will also include two other noise models discussed above: the multiplicative (6) and the truncated additive (10) one. The pure additive model (8) can be obtained from the regularized or truncated additive ones by setting the regularization or truncation parameter to zero. The regularized whitenoise model is now the default choice in PlanetPack, because so far we noticed no bad side or tradeoff consequences from its usage: the fit robustness is either improved or left unchanged. In view of that, the additive model may now be considered as deprecated. The model of the red noise in PlanetPack is now corrected to disallow negative values of the parameter in (18). PlanetPack 2.0 can also deal with the modulated red noise (24), although this is still an experimental feature. To achieve the best performance on modern multicore CPUs, some computationallyheavy PlanetPack algorithms can now be parallelized in multiple threads.
In addition to the new features that are based on the methods from the present paper, PlanetPack 2.0 will include a new optimized algorithm of evaluating the Keplerian periodogram of the RV data (Cumming, 2004), equipping it with a new fast analytic method of estimating its significance levels (Baluev, 2014). Finally, it includes a maximumlikelihood algorithm for fitting exoplanetary transit lightcurves using the new accurate limbdarkening model by Abubekerov & Gostev (2013). Although still rather experimental, this transit fitting module inherits most noise modelling features of the RV fitter, including an option of handling the red noise in the photometry data with the model (18) and the use of any of the three white noise models mentioned above.
Acknowledgements
This work was supported by the Russian Foundation for Basic Research (projects No. 120231119 mol_a and 140292615 KO_a), by the President of Russia grant for young scientists (MK733.2014.2) and by the programme of the Presidium of Russian Academy of Sciences “Nonstationary phenomena in the objects of the Universe”. I would like to express my gratitude to the reviewer, G. AngladaEscudé, for his invaluable suggestions concerning the draft of the paper.
References
 Abreu et al. (2012) Abreu J. A., Beer J., FerrizMas A., McCracken K. G., Steinhilber F., 2012, A&A, 548, A88
 Abubekerov & Gostev (2013) Abubekerov M. K., Gostev N. Y., 2013, MNRAS, 432, 2216
 AngladaEscudé & Tuomi (2012) AngladaEscudé G., Tuomi M., 2012, A&A, 548, A58
 Baluev (2008) Baluev R. V., 2008, MNRAS, 385, 1279
 Baluev (2009) Baluev R. V., 2009, MNRAS, 393, 969
 Baluev (2011) Baluev R. V., 2011, Celest. Mech. Dyn. Astron., 111, 235
 Baluev (2013a) Baluev R. V., 2013a, Astronomy & Computing, 2, 18
 Baluev (2013b) Baluev R. V., 2013b, Astronomy & Computing, 34, 50
 Baluev (2013c) Baluev R. V., 2013c, MNRAS, 429, 2052
 Baluev (2014) Baluev R. V., 2014, MNRAS, submitted, arXiv:1409.6115
 Baluev & Beaugé (2014) Baluev R. V., Beaugé C., 2014, MNRAS, 439, 673
 Butler et al. (2006) Butler R. P., et al., 2006, ApJ, 646, 505
 Cameron & Schüssler (2013) Cameron R. H., Schüssler M., 2013, A&A, 557, A83
 Chicheportiche & Bouchaud (2012) Chicheportiche R., Bouchaud J.P., 2012, Phys Rev E, 86, 041115
 Cumming (2004) Cumming A., 2004, MNRAS, 354, 1165
 Cumming et al. (1999) Cumming A., Marcy G. W., Butler R. P., 1999, ApJ, 526, 890
 Dumusque et al. (2011) Dumusque X., et al., 2011, A&A, 535, A55
 Dumusque et al. (2012) Dumusque X., et al., 2012, Nature, 491, 207
 Endl et al. (2012) Endl M., et al., 2012, ApJ, 759, 19
 Feroz & Hobson (2014) Feroz F., Hobson M. P., 2014, MNRAS, 437, 3540
 FerrazMello (1981) FerrazMello S., 1981, AJ, 86, 619
 Fischer et al. (2008) Fischer D. A., et al., 2008, ApJ, 675, 790
 Fischer et al. (2014) Fischer D. A., Marcy G. W., Spronck J. F. P., 2014, ApJS, 210, 5
 Ford (2006) Ford E. B., 2006, ApJ, 642, 505
 Gregory (2007a) Gregory P. C., 2007a, MNRAS, 374, 1321
 Gregory (2007b) Gregory P. C., 2007b, MNRAS, 381, 1607
 Hatzes (2013) Hatzes A. P., 2013, ApJ, 770, 133
 Henry et al. (2000) Henry G. W., Baliunas S. L., Donahue R. A., Fekel F. C., Soon W., 2000, ApJ, 531, 415
 Lomb (1976) Lomb N. R., 1976, Ap&SS, 39, 447
 Marcy et al. (2002) Marcy G. W., Butler R. P., Fischer D. A., Laughlin G., Vogt S. S., Henry G. W., Pourbaix D., 2002, ApJ, 581, 1375
 Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
 McArthur et al. (2004) McArthur B. E., et al., 2004, ApJ, 614, L81
 Meschiari et al. (2011) Meschiari S., Laughlin G., Vogt S. S., Butler R. P., Rivera E. J., Haghighipour N., Jalowiczor P., 2011, ApJ, 727, 117
 Naef et al. (2004) Naef D., Mayor M., Beuzit J. L., Perrier C., Queloz D., Sivan J. P., Udry S., 2004, A&A, 414, 351
 Nelson et al. (2014) Nelson B. E., Ford E. B., Wright J. T., Fischer D. A., von Braun K., Howard A. W., Payne M. J., Dindar S., 2014, MNRAS, 441, 442
 O’Toole et al. (2008) O’Toole S. J., Tinney C. G., Jones H. R. A., 2008, MNRAS, 386, 516
 Robertson et al. (2013) Robertson P., Endl M., Cochran W. D., DodsonRobinson S. E., 2013, ApJ, 764, 3
 Robertson et al. (2014) Robertson P., Mahadevan S., Endl M., Roy A., 2014, Science, 25, 440
 Saar & Donahue (1997) Saar S. H., Donahue R. A., 1997, ApJ, 485, 319
 Saar et al. (1998) Saar S. H., Butler R. P., Marcy G. W., 1998, ApJ, 498, L153
 Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
 Schneider (1995) Schneider J., 1995, The Extrasolar Planets Encyclopaedia, www.exoplanet.eu
 Tan & Cheng (2013) Tan B., Cheng Z., 2013, Ap&SS, 343, 511
 Tan et al. (2013) Tan X., Payne M. J., Lee M. H., Ford E. B., Howard A. W., Johnson J. A., Marcy G. W., Wright J. T., 2013, ApJ, 777, 101
 Tuomi & AngladaEscudé (2013) Tuomi M., AngladaEscudé G., 2013, A&A, 556, A111
 Tuomi & Jenkins (2012) Tuomi M., Jenkins J. S., 2012, arXiv.org, 1211.1280
 Tuomi et al. (2014) Tuomi M., AngladaEscude G., Jenkins J. S., Jones H. R. A., 2014, MNRAS, submitted, arXiv:1405.2016
 Wittenmyer et al. (2009) Wittenmyer R. A., Endl M., Cochran W. D., Levison H. F., Henry G. W., 2009, ApJS, 182, 97
 Wright (2005) Wright J. T., 2005, PASP, 117, 657
 Zechmeister & Kürster (2009) Zechmeister M., Kürster M., 2009, A&A, 496, 577
Appendix A A weighted analogue of the KolmogorovSmirnov test for tail probabilities
Consider a sequence of random numbers sampled from some parent distribution.^{3}^{3}3Within this appendix these designations are unrelated to the number of observations and measurements in a time series. Our task is to test whether this observed set of is consistent with some proposed cumulative distribution function (hereafter CDF) or not. This can be verified by means of the classic KolmogorovSmirnov (hereafter KS) test based on the maximum deviation between the empiric CDF , which is constructed from , and the proposed CDF :
(25) 
Under assumption that is continuous and indeed represents the CDF of , and for , the distribution function of the KS test statistic is given by the following decomposition, first obtained by A.N. Kolmogorov: