A&A 480, 563-571 (2008)
DOI: 10.1051/0004-6361:20079000

Extrasolar planet detection by binary stellar eclipse timing: evidence for a third body around CM Draconis

H. J. Deeg1 - B. Ocaña1,2 - V. P. Kozhevnikov3 - D. Charbonneau4 - F. T. O'Donovan5 - L. R. Doyle6


1 - Instituto de Astrofísica de Canarias, C. Via Lactea S/N, 38205 La Laguna, Tenerife, Spain
2 - Instituto de Radio Astronomía Milimétrica (IRAM), Av. Divina Pastora 7, Núcleo Central, 18012 Granada, Spain
3 - Astronomical Observatory, Ural State University, Lenin ave. 51, Ekaterinburg, 620083, Russia
4 - Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA
5 - California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
6 - SETI Institute, 515 N. Whisman Road, Mountain View, CA 94043, USA

Received 5 November 2007 / Accepted 28 December 2007

Abstract
Aims. Our objective is to elucidate the physical process that causes the observed observed-minus-calculated (O-C) behavior in the M4.5/M4.5 binary CM Dra and to test for any evidence of a third body around the CM Dra system.
Methods. New eclipse minimum timings of CM Dra were obtained between the years 2000 and 2007. The O-C times of the system are fitted against several functions, representing different physical origins of the timing variations.
Results. Using our observational data in conjunction with published timings going back to 1977, a clear non-linearity in O-C times is apparent. An analysis using model-selection statistics gives about equal weight to a parabolic and to a sinusoidal fitting function. Attraction from a third body, either at large distance in a quasi-constant constellation across the years of observations or from a body on a shorter orbit generating periodicities in O-C times is the most likely source of the observed O-C times. The white dwarf GJ 630.1B, a proper motion companion of CM Dra, can however be rejected as the responsible third body. Also, no further evidence of the short-periodic planet candidate described by Deeg et al. (2000, A&A, 358, L5) is found, whereas other mechanisms, such as period changes from stellar winds or Applegate's mechanism can be rejected.
Conclusions. A third body, being either a few-Jupiter-mass object with a period of 18.5 $\pm$ 4.5 years or an object in the mass range of $1.5~M_{\rm jup}$ to $0.1~M_{\odot}$ with periods of hundreds to thousands of years is the most likely origin of the observed minimum timing behavior.

Key words: stars: individual: CM Dra - stars: binaries: eclipsing - eclipses - stars: planetary systems

1 Introduction

CM Dra (LP 101.15, G225-067, GJ 630.1) is a detached spectroscopic eclipsing M4.5/M4.5 binary with one of the lowest known total masses, of $0.44~M_{\odot}$. With its nearly edge-on inclination of 89.59 $\hbox{$^\circ$ }$ (see Kozhevnikova et al. 2004 for the most recent orbital and physical elements) it was chosen as the target of the first photometric search for planetary transits. Performed by the ``TEP'' project, with an intense observing campaign during the years 1994-1999 (Doyle et al. 2000; Deeg et al. 1998), over 1000 h of coverage of that system were obtained with several 1 m-class telescopes. This lightcurve was initially searched for the presence of transits from planets in circumbinary ``P-type'' orbits with 5-60 day periods, with a negative result (Doyle et al. 2000). The same lightcurve provided, however, a further possibility to detect the presence of third bodies, from their possible light-time effects on the binary's eclipse minimum times.

To date, no unambiguously circumbinary planets have been detected[*], and their discovery would constitute a new class of planets. Motivation for this work also arises from previous successes of precise timing measurements to detect the presence of planets. The first known extrasolar planets were detected through light-time effects in the signals of the Pulsar PSR 1257+12 (Wolszczan & Frail 1992) and recently, sinusoidal residuals in the pulsation frequency of the sdB pulsating star V391 Peg have been explained through the presence of a giant planet (Silvotti et al. 2007), leading to the first detection of a planet orbiting a post-red-Giant star. In both cases, timing measurements have led to detections of planets that would have been difficult or impossible to find with other planet-detection methods, a situation that is similar to the detection of planets around eclipsing binaries - unless they exhibit transits.

A first analysis of 41 eclipse minima times for CM Dra, presented in Deeg et al. (2000) (hereafter DDK00) gave a low-confidence indication for the presence of a planet of 1.5-3 Jupiter masses, with a period of 750-1050 days. This result was based on a power-spectral analysis of the minimum timings' residuals against a linear ephemeris, indicating a periodic signal with an amplitude of about 3 s. Motivated by the result of DDK00, we continued surveying the CM Dra system with occasional eclipse observations in the following years. During this time, it became increasingly clear that a simple linear ephemeris would not provide a sufficient description of the general trend of the eclipse times any longer. This led to the objective of this paper - a thorough and systematic discussion of the possible processes acting on this interesting system, and an evaluation of the presence of a third body of planetary mass or heavier.

In the following sections, we present the data for this analysis (Sect. 2), evaluate the effect of several physical processes on the eclipse timings (Sect. 3), and compare the significance of several numerical fits for the explanation of the observed trend (Sect. 4). This is followed in Sect. 5 by a discussion of the physical implications of these fits, with conclusions in Sect. 6.

  
2 The observational data

The minimum times analyzed here include all minimum timings derived in DDK00. Since its publication, we performed dedicated eclipse observations with the IAC80 (0.8 m) and the INT 2.5 m telescopes within the Canary Islands' Observatories and with the Kourovka 0.7 m telescope of the Ural State University in Ekaterinburg, Russia. In summer 2004, one of the fields surveyed for planetary transits by the TrES collaboration (Alonso et al. 2007) contained CM Dra and time series spanning several weeks were obtained from the Sleuth 10 cm telescope (O'Donovan et al. 2004). In all cases, photometric time series were obtained from which the eclipse times were measured using the method of Kwee & van Woerden (1956). We accepted only results where the formal measurements error from that algorithm was less than 10 s, which required data with a largely complete and uninterrupted coverage of individual eclipses. The timings obtained were then corrected to solar-system barycentric (BJD) times with the BARYCEN routine in the ``aitlib'' IDL library[*] of the University of Tübingen.

Only a few timings of CM Dra of comparable quality could be found in the literature: Lacy (1977) gives a total of 9 minimum times from observations in 1976. These timings scatter in O-C by about 30 s, and were based on photoelectric data with frequent gaps, including several incomplete eclipses. We therefore chose to remeasure them from Lacy's tabulated lightcurve with the same procedures used for our own data, from which only two timings of sufficient quality could be obtained for the further analysis. Only two additional timings from the years 2004 and 2005 could be found in the literature (Dvorak 2005; Smith & Caton 2007), both in good agreement with the other data. The final data set contains 63 minima, of which 27 are primary and 36 are secondary eclipses.

Since these minimum timings cover about 30 years of observations and are of consistently high precision, the effects from the 18 leap-seconds that have been introduced into Universal Time during that span need to be corrected for. Consequently, all minimum times used in this work were converted from the conventional UT to TAI (International Atomic Time)[*] which is a timescale with constant and uniform flow, without discontinuities from leap-seconds. All these minimum times are listed in Table 1 together with the formal errors of the minimum times, an indication (I/II) for primary or secondary eclipses, the cycle number E and O-C (observed - calculated) residuals against the ephemerides of Deeg et al. (1998). For data newly presented in this work, the originating telescope is also indicated (``Kourv.'' refers to the Kourovka telescope).

Table 1: CM Dra eclipse minimum times, given in International Atomic Time (TAI).

Table 2: Parameters of fits to O-C times.


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{9000fig1.eps} \end{figure} Figure 1: O-C values of CM Dra eclipse minimum times from Table 1, against ephemerides from DDK00.
Open with DEXTER

  
3 Physical processes' effects on O-C times

The further analysis of the eclipse minimum times is based on the analysis of the temporal development of the ``O-C'' residuals between observed and calculated (expected from ephemeris) eclipse minimum times. In these,

\begin{displaymath}% ({\rm O{-}C})_{E} = T_{E} - T_{{\rm c},E}, \end{displaymath} (1)

where $T_{{\rm c},E}$ refers to a minimum time calculated from an ephemerides at ``Epoch'' or cycle number E, and TE refers to the corresponding observed minimum time. The observed times TE are related to the minimum times T'E in the binary system's rest frame by TE = T'E + dE/c, where dE is the distance from the observer to the binary at cycle E and c is the velocity of light. Hence, based on a linear ephemeris $T_{{\rm c},E} = E P_{\rm c} + T_{\rm c,0}$, where $P_{\rm c}$ and $T_{\rm c,0}$ are the ephemerides' period and time of conjunction. The difference, $({\rm O{-}C})_{E}$, is given by:
 
$\displaystyle % ({\rm O{-}C})_{E} = T'_{E} + d_{E}/c - E P_{\rm c} - T_{\rm c,0}.$     (2)

This general expression allows us now to develop the cases to be considered in our analysis.

  
3.1 Accelerated binary systems with constant period

First, we review the case of a binary system moving at a constant velocity v0 relative to the observer, with $v_0 \ll c$. The distance to this system is given by dE = d0 + v0 E P', and the minimum times in the binary system's rest-frame are given by T'E = T'0 + E P', where P' is constant. The sub-indices ``0'' refer to the parameters' values at the moment when E=0. From Eq. (2), we obtain then a relation linear in E:

 
                      $\displaystyle % ({\rm O{-}C})_{E}$ = $\displaystyle T'_0 + E P' + (d_0 + v_0 E P')/c - E P_{\rm c} - T_{\rm c,0}$ (3)
  = $\displaystyle (T_0 - T_{\rm c,0}) + E(P - P_{\rm c})$ (4)

where P = P'(1+v0/c) is the observable period and T0 = T'0 + d0/c is the minimum time that should have been observed at E=0. Hence, in a system where the observed O-C times deviate linearly from the ephemeris, the terms $\kappa_0= (T_0 - T_{\rm c,0})$ and $\kappa_1=(P - P_{\rm c})$ indicate errors in the derivation of the original ephemerides $(T_{\rm c,0},P_{\rm c})$, but do not have any further physical meaning.

If we consider as D(E) any additional distance to the binary that cannot be expressed by a linear term, so that dE = d0 + v0 E P' + D(E), then the O-C times will be modified by the ``light-time effect'' D(E)/c:

 
$\displaystyle % ({\rm O{-}C})_{E} = (T_0 - T_{\rm c,0}) + E (P - P_{\rm c}) + D(E)/c.$     (5)

Hence, the non-linear components of $({\rm O{-}C})_{E}$ describe the non-linear components of the distance to the binary. The non-linear distance component D has necessarily to be caused by an acceleration process, with $D(E) = \int\int_{T'_0}^{T'_0 + E P'} a_\parallel(t)~ {\rm d}t {\rm d}t$, where $a_\parallel$ is the acceleration along the line of sight to the binary. Within the scope of this work, the following acceleration scenarios have been considered:
i)
For the case of a constant acceleration  $a_\parallel$, with $D(E) =\frac{1}{2}a_\parallel (EP')^2$ we obtain now:
 
                   $\displaystyle % ({\rm O{-}C})_{E}$ = $\displaystyle (T_0 - T_{\rm c,0}) + E (P - P_{\rm c}) + E^2 a_\parallel \frac{P'^2}{2c}$ (6)
  = $\displaystyle \kappa_0 + E\kappa_1 + E^2\kappa_2$ (7)

where the $\kappa_i$ refer to the polynomial coefficients of the fit given in Table 2.

ii)
In the case of an acceleration that undergoes a constant change $\dot{a}_\parallel$, the O-C times behave like a third order polynomial:
$\displaystyle % ({\rm O{-}C})_{E} = \kappa_0 + E\kappa_1 + E^2\kappa_2 + E^3\kappa_3$     (8)

where $\kappa_3$ is given by:
$\displaystyle % \kappa_3= \dot{a}_\parallel \frac{P'^3}{6c}$     (9)

whereas the coefficients $\kappa_0 - \kappa_2$ remain identical to the previous case.

iii)
For a binary that is accelerated due to third body on a circular orbit, the amplitude of the timing variation D/c is then given by:

 \begin{displaymath}% D/c = \frac{m_3\ d_{\parallel {\rm b3}}}{(m_{\rm b} + m_3)\... ...3)^{2/3}c} \left( \frac {P_3^2 G}{4\pi^2} \right)^{1/3} \sin i \end{displaymath} (10)

where $d_{\parallel {\rm b3}}$ is the line-of-sight-component of the distance between the barycenter of the binary and the third body, and m3 and $m_{\rm b}$ are the masses of the third body and the binary stars, respectively. In the right hand term, $d_{\parallel {\rm b3}}$ has been substituted for P3, the period of the third body, with i being the inclination of its orbital plane and G is the gravitational constant. A development of the general case is given by Irwin (1959), with further examples of recent applications in Demircan & Budding (2003).

From Eq. (10), the O-C times are then given by

 
$\displaystyle % ({\rm O{-}C})_{E} = \kappa_0 + E\kappa_1 + \kappa_{\rm d} \cos ( E \kappa_\phi+ \phi_0)$     (11)

with $\kappa_{\rm d} =D/c$, and the term $E \kappa_\phi+ \phi_0$ describing the phase of the third body, where

 \begin{displaymath}% \kappa_\phi = \frac{P' 2 \pi}{P_3} \end{displaymath} (12)

and $\phi_0$ is the phase of the third body at T0.

  
3.2 Variation of the intrinsic binary period

Here we consider the consequences of a period variation of a binary moving at a constant velocity. The intrinsic period variation  $\partial P'/\partial E$ shall be small enough to be considered constant across the observed time span. The period at cycle E is then given by $P'_{E} = \int \frac{\partial P'}{\partial E} ~{\rm d}E + P' _{0} = E \frac{\partial P'}{\partial E} + P' _{0}$, and the times of minima in the binary's reference frame are:

                  T'E = $\displaystyle \int P'_{E}~{\rm d}E + T' _{0}$ (13)
  = $\displaystyle \frac{1}{2}E^2 \frac{\partial P'}{\partial E} + E P' _{0} + T' _{0}.$ (14)

After converting from the times T'E to the observable times TE by inserting this equation into Eq. (2), we obtain a quadratic equation that is similar to the accelerated system described by Eq. (7), with the only difference being that the parameter $\kappa_2$ is now given by:

 \begin{displaymath}% \kappa_2= \frac{1}{2} \frac{\partial P'}{\partial E}\left(1+\frac{v_0}{c}\right)\cdot \end{displaymath} (15)

   
4 The observed minimum times

As in DDK00, O-C times were derived using the linear ephemerides given by Deeg et al. (1998), which was based on a fit to eclipse timings observed from 1994 to 1996. For the present work, this ephemerides was converted to TAI by adding 29 s, which corresponds to the difference TAI-UT that was in effect during most of that time (July 1, 1994 to Dec. 31, 1995). The O-C values of DDK00 differ therefore by a maximum of only 1 s between the original work in DDK00 and the present work. The ephemerides conversion from Deeg et al. (1998) to TAI is:

    $\displaystyle T_{{\rm cI}} = 2~449~830.75734 \pm 0.000~01 + P_{\rm c} \ E\;({\rm TAI})$ (16)
    $\displaystyle T_{{\rm cII}} = 2~449~831.39037 \pm 0.000~01 + P_{\rm c} \ E\;({\rm TAI})$ (17)

where $T_{{\rm cI}}$ and $T_{{\rm cII}}$ refer to primary and secondary minima, respectively, E is the epoch or cycle number, and the period is given by $P_{\rm c} = 1.268~389~861$ $\pm$ 0.000 000 005 days. The O-C values against that ephemerides are given in Table 1 and shown in Fig. 1. Since the writing of DDK00, a clear trend of increasing O-C values has become apparent, which is also apparent from the extension into the past through the inclusion of the values from Lacy (1977), which weren't included in DDK00's analysis. The linear ephemerides of DDK00, therefore, no longer provides the best description of the observed O-C values. However, we chose to maintain this ephemerides, since small errors in the parameters of a linear ephemerides have no physical meaning in the interpretation of higher-order O-C dependencies, as was shown in Eq. (4).

4.1 Temporal evolution of the phase of the secondary eclipse

A primary concern was assuring that the primary and secondary eclipses were not undergoing any evolution of their relative phase due to, for example, variations in eccentricity or in the argument of periastron due to apsidal motion. Therefore we investigated if the orbital phase of the secondary eclipse underwent any variations. This was done through a comparison of two sub-samples of the timings from Table 1, taking an early one from Epochs 0 to 400 and a late one from Epochs 2500-2900. In both samples, the average of the primary eclipses was set to zero, and the corresponding phase of the secondary eclipses were calculated, giving:

\begin{displaymath}% {\rm phase }=0.499064 \pm 0.000017 {\rm\ for}\ E=0{-}400 \end{displaymath} (18)


\begin{displaymath}% {\rm phase }=0.499039 \pm 0.000029 {\rm\ for}\ E= 2500{-}2900. \end{displaymath} (19)

The difference between these two values is well within their error-bars; hence no relevant change in the phase of the secondary eclipse was detected. Since primary and secondary eclipses of CM Dra have very similar depths, resulting in timing measurements of similar precision, both types of eclipses were treated together and equally in the further analysis.

  
4.2 Model fits to the O-C times

Our initial intent was a direct fitting of the functions described in Sect. 3 to the O-C residuals of Table 1. However, in the subsequent statistical analysis we noted two problems: First, the average formal error in the individual O-C measurements from Table 1 is 3.90 s. On the other hand, the standard deviation among the residuals (observations - fits) could not be reduced below 5.6 s. Reduced $\chi ^2$ values were not less than $\chi_{\rm red}^2 \approx 4.2$, even when applying 5th or 6th order polynomial fits, whereas a ``good'' fit should indicate values of $\chi_{\rm red}^2 \approx 1$. Since the fits are not intended to - and cannot - model the point-to-point variations and furthermore, since the presence of periodic short-frequent O-C timing variations can be ruled out (see Sect. 5.3), we concluded that the formal errors are sub-estimating the real measurement errors. In the further analysis, a value of 5.6 s was adopted as a minimum measurement error and errors smaller than this were set to 5.6 s. With fourth to sixth order polynomial fits to these data we now obtain $\chi_{\rm red}^2 \approx 1$.

A second problem arose because the data-points aren't uniformly or randomly distributed, but clustered into yearly observing seasons. These clusters act like pivots for the fitting functions, reducing the degrees of freedom of the model fits over the number of data points. The choice of the correct number of degrees of freedom is, however, important for the model comparison performed in Sect. 4.4. Consequently, for each year of observations, we generated one single data point from a weighted average of each season's points. The resulting binned O-C times are shown in Fig. 2.


  \begin{figure} \par\includegraphics[width=7.6cm,clip]{9000fig2.eps} \end{figure} Figure 2: Seasonally averaged O-C values of CM Dra, with parabola fit (dashed line) and sine-linear fit (solid line).
Open with DEXTER

Fits of the functions that have been discussed in Sect. 3 were then performed against this sample, with their best-fit parameters given in Table 2. The fits were obtained with the IDL library routines ``POLY_FIT'' for the polynomial fits, and ``AMOEBA'' (based on Press et al. 1992) for the sinusoidal fit, with both algorithms performing minimizations of the $\chi ^2$ values. It should be noted that the aim of the sine-linear fit was to test if the data can be well fit to one or a few cycles of such a function, with a period longer than about 10 years. For a search for shorter periodicities see Sect. 5.3.

The parameters' errors given in Table 2 are 1-sigma confidence limits, whose derivation is described in the next section. Furthermore, the table's last three columns give the best-fit $\chi ^2$ values, calculated for a sample standard deviation of s=3.2 s, the Akaike Information Coefficient AIC$_{\rm c}$, and the Akaike weights wi; both of which will be introduced in Sect. 4.4.

4.3 Errors of fit-parameters

For an estimation of the errors of the fit-parameters we applied initially the common resampling or bootstrap method (e.g. Cameron et al. 2007), consisting of repeated model fits to a synthetic data set. It requires one to assume some reference function that describes correctly the general trend of the data, against which residuals of the data points are generated. The synthetic data are generated by permutating the residuals among the data points. The model to be evaluated is then fitted against the synthetic data and distributions of the obtained fit-parameters are used to estimate the likelihood distribution of the parameters from the model-fit on the original data. This method is easy to implement and leads to synthetic samples with properties similar to the original data, but in this work's context two problems arose:

First, the analysis is based on the assumption that the distribution of the fit-parameters on the synthetic sets is identical to the likelihood distribution of the parameters obtained from the fit of the original data, which is far from certain.

Second, resampling is based on the assumption that the reference function is a correct description of the trend of the data, and that residuals against it are measurement errors. This approach may be justified if the underlying physical model - and the function that describes it - is known, and only a refining of parameters is required. In this work however, we are also faced with an uncertainty about the nature of the reference function.

A method that overcomes these problems, and which has come to the awareness of astrophysicists in recent years, is the Markov Chain Monte-Carlo (MCMC) method (e.g. Holman et al. 2006; Burke et al. 2007; Tegmark et al. 2004; Ford 2005). MCMC is based on the states of a Markov Chain undergoing random variations, whose probabilities are however directed by the maximum likelihood estimator (MLE) of the fit-parameters (relative to the original data) at each step in the chain. The frequency of states of the Markov Chain at a given parameter value indicates then the posterior probability distribution of that parameter. We refer to Ford (2005) for further references to the MCMC, as well as for the implementation of the MCMC with the Metropolis-Hastings Algorithm that was used in our data analysis.


  \begin{figure} \par\includegraphics[width=8cm,clip]{9000fig3.eps} \end{figure} Figure 3: Distribution of the $\kappa _{\phi }$ parameter from the sine-linear model against $\chi ^2$ from a Markov Chain of 20 million steps. The dashed lines indicate the parameter's 1-$\sigma $ confidence region, corresponding to $\Delta \chi ^2 = 1$ over the best fit's $\chi ^2$ value (cross).
Open with DEXTER

Our final error-analysis is however not based on the posterior probability distributions derived from histograms of the parameter distributions. Instead we have used the MCMC as a tool to explore the relation of $\chi ^2$ against the multivariate parameters. Similar to Burke et al. (2007), for any recorded step of the MCMC chain, the encountered parameters and the corresponding $\chi ^2$ values were registered. As an example, Fig. 3 shows the relation between the $\kappa _{\phi }$ parameter in the sine-linear model fit against $\chi ^2$. There is a clear lower limit of $\chi ^2$ for a given parameter value. With increasing numbers of iterations in the Markov Chain, this lower limit approaches the best possible fit at a given parameter value. Hence, the lower limits of $\chi ^2$ may be used as tracings of the best-fit $\chi ^2$ against one or multiple parameters. Since the distribution-density of a MCMC sequence gravitates towards the regions of lowest $\chi ^2$, the MCMC method may therefore be used as a simple tool to trace low-sigma confidence regions around the best-fit parameters. This application of the MCMC also avoids a problem that easily occurs in the interpretation of posterior probabilities from the parameter densities. As can be seen in Fig. 3, the distribution of points also has dense zones close to the left cutoff, although the minimum value of $\chi ^2$ in that zone corresponds to very poor fits. With the flat dependency of the best-fit $\chi ^2$ against $\kappa _{\phi }$ in that zone, the Markov Chains had difficulties returning to better-fitting values. In the given case, the Markov Chains instead went into ``exploring'' a wide multi-parameter space of poor models that opened up close to the left cutoff in Fig. 3.

For the errors indicated in Table 2 we used chains with a length of 20 million iterations, after discarding the first 1000 steps of the chains.

  
4.4 Model selection

After applying the fits of several models, physical interpretations will need some information about the likelihood that a given fit corresponds to the true behavior of the observed data. This question is commonly referred to as ``model selection''. Values such as $\chi ^2$ (see Table 2), or the ``reduced Chi-square'' of $\chi^2 / \nu$, where $\nu$ is the degrees of freedoms, give a general indication of the quality of a fit. However, in the case of small or moderate differences of fit-quality among models, they serve little to generate statements that are useful for the model selection.

Probably the best-known method of assigning likelihoods in fit-comparisons is the ``F-test''. It is however valid only for the comparison of nested models, such as polynomials of different orders, and was therefore not used further. The Akaike Information Coefficient (AIC, Akaike 1974) does not have this limitation, which led to its use in the further investigation. For an introduction to the use of the AIC as a tool for model selection, we refer the reader to Burnham & Anderson (2004); Liddle (2007); Mazerolle (2004). Using the residuals' squared sum (RSS) as the likelihood estimator, where ${\rm RSS} = \sum{(y_i - f_i)^2}$ with yi being the data and fi the model values, the AIC is given by:

\begin{displaymath}% {\rm AIC} = -n \ln~ ({\rm RSS}/n) + 2 k \end{displaymath} (20)

where n is the number of observations and k is the number of model parameters. Here we use the generally preferred (Burnham & Anderson 2004) second order corrected coefficient ``AIC$_{\rm c}$'', which is valid for both small and large samples:

\begin{displaymath}% {\rm AIC}_{\rm c} = {\rm AIC} + \frac{2k(k+1)}{N-k-1}\cdot \end{displaymath} (21)

While smaller AIC values indicate better fits, their absolute values don't have any meaning. They are only useful if values from different models are compared. If the best among several models has a value of ${\rm AIC}_{{\rm c,~min}}$, then for any model i the differences $\Delta_i = {\rm AIC}_{{\rm c},i} - {\rm AIC}_{{\rm c,~min}}$ may be calculated. Differences of $\Delta_i \le 2$ indicate substantial support (evidence) for model i; models where $4 \le \Delta_i \le 7$ have considerably less support, while models with $\Delta_i \ge 10$ have essentially no support (Burnham & Anderson 2004). For the comparison within a set of R models, normalized Akaike weights may be derived for each model i, with

\begin{displaymath}% w_i = \frac{\exp~(-\Delta_i/2)}{\sum_{r=1}^R \exp~ (-\Delta_r/2)} \end{displaymath} (22)

where all weights wi sum up to 1 (see Table 2). Following Akaike (1981), these weights may be interpreted as a likelihood that can be assigned to each of the models, with the parabolic fit being most likely, followed closely by the sine-linear fit. A word of caution, also reflected in several references about this topic (e.g. Liddle 2007), should however be given against the use of these statistical values as a strong argument in favor of one or the other model: there are several alternative indicators available, such as the Bayesian Information criterium (BIC, Schwarz 1978) or the Deviance Information Criterion (DIC, Spiegelhalter et al. 2002), with different ``penalizations'' for models with additional degrees of freedom. While these criteria indicate similar preferences for models with well separated $\chi ^2$ or AIC$_{\rm c}$ values, interchanges in ranking may happen among models that are close in AIC$_{\rm c}$ values. This cautionary position was backed by a calculation of the BIC for our models. In that case, the sine-linear fit was ranked best, with a slightly lower (and hence ``better'') BIC than the parabola fit, while the ranking of the other models remaining unchanged.

In summary, both the simple linear fit and the third order polynomial fit have significantly less support than the parabolic fit and the sine-linear fit. In Sect. 5 we therefore focus on the physical implications from these two top-ranked models.

  
5 Discussion: possible causes of the observed O-C times

Common to all fits, the linear parameter $\kappa_1$ has fairly similar values indicating clearly that the average period of CM Dra across the recorded observations is several milliseconds longer than given by Deeg et al. (1998). As shown in Sects. 3.1 and 3.2, a parabolic O-C function may arise from two causes, an intrinsic variation in the system's period, or a light-time effect from a constant acceleration of the entire binary system. In both cases, the only interesting parameter is the quadratic term, found to be (see Table 2) $\kappa_2 = (7.0 \pm 1.3)$ $\times $ 10-7 s/period.

5.1 Intrinsic period variation

Considering an intrinsic period variation, the change in period-length per cycle is given by Eq. (15) with $\frac{v_0}{c_l} \ll 1$ as:

\begin{displaymath}% \frac{\partial P'}{\partial E} \approx 2 \kappa_2 = (1.4\pm0.26)\times10^{-6}~{\rm s/period}, \end{displaymath} (23)

The corresponding unitless period change per time is given by:

 \begin{displaymath}% \frac{\partial P'}{\partial t}=\frac{\partial P'}{\partial E} \frac{1}{P'}=1.28\times10^{-11}. \end{displaymath} (24)

Demircan et al. (2006) performed a statistical study of the orbital parameters of a sample of detached chromospherically active eclipsing binaries of different ages. Out of that sample, they concluded that their periods decrease with an average value of $\alpha=3.96$ $\times $ $10^{-10}~{\rm yr}^{-1}$, with $\alpha$ defined by the differential equation ${\rm d}P/{\rm d}t = - \alpha P$. This value may be considered constant throughout a large part of a binary's evolution and is due to angular momentum loss from magnetically driven stellar winds. The corresponding period change of CM Dra would be $\frac{\partial P'}{\partial t}=-1.38$ $\times $ 10-12, obtained by multiplying $-\alpha$ with CM Dra's period in units of years (3.47 $\times $ 10-3 yr). This value is of opposite sign than the observed one (Eq. (24)), and is an order of magnitude smaller. Hence, angular momentum loss may well be present, but is not detectable in the current minima timings.

5.2 Acceleration due to a quasi-stationary attractor

The second source for a parabolic shape in an O-C diagram could be a constant acceleration of the binary along the line of sight, with changes of acceleration strength and direction during the span of observations being negligible. The acceleration term is given from Eq. (7) as:

 \begin{displaymath}% a_\parallel = 2 \frac{c \kappa_2}{P^{'2}} = (3.5\pm0.7)\times10^{-8}~{\rm m/s}^2 \end{displaymath} (25)

or

\begin{displaymath}% \frac{a_\parallel}{c}= (1.17\pm0.22)\times10^{-17}~{\rm s}^{-1}. \end{displaymath} (26)

We note that this acceleration is at least two orders of magnitude larger than the acceleration of the Solar system, currently constrained within a few $\times 10^{-19}~{\rm s}^{-1}$ (Zakamska & Tremaine 2005). A quasi-constant acceleration may be caused by a third body of mass m3 at a distance far away enough so that mutual orbital motions don't lead to significant changes in the acceleration vector. The acceleration on the binary caused by a third body is given by:

 \begin{displaymath}% \vec{a}= \frac{Gm_3}{r^3} \vec{r} \end{displaymath} (27)

where $\vec{r}$ is a distance vector from the barycenter of the binary towards the third body. We note that this acceleration is independent of the mass of the binary. The acceleration component along the line of sight is then given by

\begin{displaymath}% a_\parallel = \frac{Gm_3}{r_\perp^2}\cos^2 i \sin i \end{displaymath} (28)

where i is the inclination, defined here as the angle between $\vec{r}$ and the plane of the sky, and $r_\perp= r \cos i$ is the lateral component of r. The equation above gives $a_\parallel = 0$ for inclinations of both $0\hbox{$^\circ$ }$ and $90\hbox{$^\circ$ }$ and a maximum for $a_\parallel$ at inclinations of $i = \arctan \sqrt{1/2} = 35.26\hbox{$^\circ$ }$, leading to

\begin{displaymath}% a_{\parallel} \le 0.3849\ \frac{Gm_3}{r_\perp^2}\cdot \end{displaymath} (29)

Converting to common astronomical units, the minimum mass m3 at a given lateral distance is then obtained by

 \begin{displaymath}% \left(\frac{m_3}{M_{\odot}}\right) \ge 438.26 \left(\frac{a... ...m s}^{-2}}\right) \left(\frac{r_\perp}{{\rm AU}}\right)^2\cdot \end{displaymath} (30)

Replacing $r_\perp$ by an angular separation based on the distance to CM Dra (15.93 pc; Chabrier & Baraffe 1995), and using $a_\parallel =3.5\pm0.7\times10^{-8}~{\rm m/s}^2$, we may now derive the minimum mass of a possible third body at a given angular separation from CM Dra:

 \begin{displaymath}% \left(\frac{m_3}{M_{\odot}}\right) \ge 0.0030 \left(\frac{\alpha}{{\rm arcsec}}\right)^2\cdot \end{displaymath} (31)

The third order polynomial fit, while resulting in a lower weight in the model selection, is not to be discarded completely. It would describe a system undergoing a constant change in acceleration  $\dot{a}_\parallel$. There is however no physical process that generates a truly constantly varying acceleration. Hence the third-order polynomial fit may only describe cases where $\dot{a}_\parallel$ is quasi-constant across the observing time-span. The third order polynomial fit may, however, provide the first terms of a Taylor-expansion of the true acceleration process of the system. For a slowly changing acceleration, like a cyclic one with long periods of O(100 yr) or longer, the 3rd order polynomial may therefore give a better description of the O-C times than the parabolic fit does. From our fit, however, with the value of $\kappa_2$ being very close to the one from the parabolic fit, the derived acceleration  $a_\parallel$ and the constraint for a third mass given in Eq. (31) do not significantly differ.

A possible source for an acceleration of the CM Dra system may be the nearby white dwarf GJ 630.1B (WD 1633+57, LP 101.16, G225-068). This has long been recognized as a proper-motion companion to CM Dra (Giclas et al. 1971) at an angular distance of 26 $\hbox{$^{\prime\prime}$ }$, which corresponds to a lateral distance $r_\perp$ of 414 AU. With a period of O(104) yr, the criteria of a quasi-constant acceleration vector during the 31 years of observational coverage is clearly given. Following Eq. (31), a minimum mass of m3 of 2.0 $\pm$ $0.4~M_{\odot}$ is however obtained for the white dwarf - much above the typical white dwarf masses of 0.5-0.7 $M_{\odot }$, and clearly above the Chandrasekar limit for white dwarfs of $\approx$ $1.4~M_{\odot}$. Assuming a mass of $0.6~M_{\odot}$ for GJ 630.1B, this object would contribute an acceleration of only $a \la 1$ $\times $ $10^{-8}~{\rm m/s}^2$ on the CM Dra system.

While GJ 630.1B can be ruled out as a source of the observed O-C variations, they may be caused by still undiscovered bodies in the brown-dwarf mass regime (13 to 80  $M_{\rm Jup}$) at maximum distances of about 5 arcsec, or by a planetary-mass object (with less then 13 Jupiter masses) at a maximum distance of 2 arcsec.

  
5.3 An orbiting third body

As shown in Sect. 4.4, the sinusoidal fit matches the observed O-C times about as well as the parabolic one. From the fitted parameter, $\kappa _{\phi }$, we obtain with Eq. (12) and $m_3 \ll m_{\rm b}$ a period of

 \begin{displaymath}% P_3=\frac{P' 2 \pi}{\kappa_\phi}=18.5\pm4.5~{\rm yr} \end{displaymath} (32)

Eq. (10), with $D/c = \kappa_{\rm d}$ can be rewritten as:
$\displaystyle % m_3 \sin i = \kappa_{\rm d} \left( \frac{m_{\rm b}}{M_{\odot}} \left/ \frac{P_3}{{\rm yr}}\right. \right)^{2/3} 2.1~M_{{\rm Jup}}.$     (33)

With the above value for P3 and the fitted one for $\kappa_{\rm d}$, this leads to:

\begin{displaymath}% m_3 \sin i = 1.5\pm0.5~M_{{\rm Jup}}. \end{displaymath} (34)

We note that such an object would have an orbital half-axis of about 5.3 AU, with a maximum separation from CM Dra of about 0.35 arcsec.

While the sinusoidal fit indicates a periodicity on a time-scale of 20 years, we also performed a search for higher frequencies. For this, the same sine-wave fitting algorithm used in DDK00 was employed. This analysis was performed on the O-C residuals using the parabola fit, and included only the relatively dense surveying that started in 1994, thereby excluding the two isolated early values obtained from Lacy (1977). The resulting power-spectrum is shown in Fig. 4.

The single broad peak at a period of 950 days and with an amplitude of 3 s that was apparent in Deeg et al. (1998)'s Fig. 2c has now disappeared, and been replaced by several peaks with amplitudes close to 3.5 s. We note that this amplitude is close to the sample standard deviation of the yearly averaged data of s=3.2 s. Since no single peak is outstanding, no indications for any periodicites on timescales of $\la$10 years remain.


  \begin{figure} \par\includegraphics[width=7.2cm,clip]{9000fig4.eps} \end{figure} Figure 4: Power-spectrum obtained from residuals of O-C values against the parabola fit.
Open with DEXTER

5.4 Applegate's mechanism

We also evaluated the mechanism introduced by Applegate & Patterson (1987) and Applegate (1992), which may give rise to orbital period modulations of binaries from the periodic variation of the shape - and of the quadrupole moment - of a magnetically active binary component across its activity cycles. The notion that a component of CM Dra may be magnetically active cannot be completely discarded, since several large flare events have been reported in the literature (Kim et al. 1997; Deeg et al. 1998; Lacy 1977).

In Applegate's model, the angular momentum of the entire system remains constant, but the distribution of angular momentum between the components' mutual orbit and the active star's internal rotation varies. Energy taken up by the variation in the internal rotation has to be reflected in the active star through a corresponding luminosity variation. Applegate's original calculation considers the energy taken up by differential rotation between an inner stellar core and an outer shell of 0.1 $M_{\odot }$, something which is inappropriate for either component of CM Dra, with masses of 0.207 and 0.237 $M_{\odot }$ (Kozhevnikova et al. 2004), respectively. We followed therefore the more general calculation introduced by Brinkworth et al. (2006), which leads to the rotational energy that has to be provided in order to explain a given period change, while doing this for any distributions of the stellar mass into core and shell. The period variation under consideration is given by (Applegate 1992):

\begin{displaymath}% \Delta P = P\ 2\pi \frac{\kappa_{\rm d}}{P_3} \end{displaymath} (35)

where $\kappa_{\rm d}$ is the amplitude of the O-C variation, P the binary period and P3 the modulation period. With corresponding values taken from Table 2 and Eq. (32), we obtain $\Delta P = 0.010$ s. The minimum energy to produce such a period change for any distribution of stellar mass into core and shell (assuming a mass distribution following the Lane-Emden equation for a polytrope of n=1.5) amounts to 3.8 $\times $ 1042 erg if CM Dra A is considered as the active star (see Fig. 5); it would be slightly higher (4.4 $\times $ 1042 erg) in the case of CM Dra B. This energy may be compared to the luminosity of CM Dra of 1.06 $\times $ $10^{-2}~L_{\odot}$ (Kozhevnikova et al. 2004), which corresponds to a release of radiant energy of 2.4 $\times $ 1040 erg over the same span of P3=18.5 years. With the energy required for the period change being two orders of magnitude larger than the radiant energy, Applegate's mechanism can definitively be discarded as a source of the period variations.


  \begin{figure} \par\includegraphics[width=7.15cm,clip]{9000fig5.eps} \end{figure} Figure 5: Energy required to vary the internal rotation of CM Dra A in order to reproduce the observed period change with Applegate's model, for all possible values of CM Dra A's shell mass. The lowest amount of energy required is 3.8 $\times $ 1042 erg at a shell-mass of 0.016 $M_{\odot }$.
Open with DEXTER

  
6 Conclusions

The O-C timing shows a clear indication of a non-linear trend. After a review of potential causes (previous section), the remaining explanations are given by the presence of a third body. The nearly equal statistical weight of the parabolic and the sinusoidal fits currently prevents setting a clear preference for the one or the other model. For the period of the third body there exist two distinct possibilities: a Jupiter-type planet with $M \sin i$ of 1-2  $M_{\rm jup}$ with a period of 18.5 $\pm$ 4.5 yr, or an object such as a giant planet or heavier, with a period of hundreds to thousands of years. Intermediate-length periods of about 25-100 years are less likely due to poor compatibility with either the sinusoidal or the polynomial fits. Regarding shorter periodic bodies, the power-spectrum shown in Sect. 5.3 excludes O-C modulations with amplitudes of larger than $\approx$3.5 s with periods of $\la$10 years. The sensitivity of O-C timing detections against orbiting third bodies decreases with their period, and hence Jupiter-mass objects on such shorter periods cannot be excluded. A low-significance candidate for such an object with a period of about 900 days was presented in DDK00. That candidate was based on a single peak in a power-spectrum with an amplitude of 3 s, whereas the newer data show several peaks of amplitudes up to 3.5 s. The newer data therefore cannot refute that candidate, however it has become most likely to have been the result of a fortuitous combination of O-C timing values.

For long-period bodies at a quasi-constant distance and position during the observed time-span, Eq. (31) allows us to set a maximum lateral separation from CM Dra for any given third-body mass. We also assume that any nearby body larger than about $\approx$ $0.1~M_{\odot}$ would have become apparent in existing images of the CM Dra field, for which a large collection of CCD images exist from the TEP project (Doyle et al. 2000; Deeg et al. 1998), or in 2Mass images in the IR. A maximum lateral distance of 6 arcsec, or 95 AU may therefore be set for the presence of an as yet undiscovered third body of $\la$ $0.1~M_{\odot}$. The corresponding maximum distances for undiscovered brown dwarfs or planets that could explain CM Dra's O-C timing behavior are 5 and 2 arcsec, respectively. While the setting of third-body minimum masses (and implied brightnesses) for a given lateral distance aids in the definition and interpretation of observing projects, we note however that the observed acceleration term may be caused by relatively small third-body masses. If we take 100 years as the shortest circular period that may mimic the quasi-stationary case, such an object's orbital distance would be 16.4 AU. If it is aligned such that $\vert\vec{r}\vert \approx r_{\parallel}$ (e.g. i close to $90\deg$), then solving Eq. (27) indicates a mass of about $1.5~M_{\rm jup}$. This mass may be considered the absolute minimum mass for a third body at a quasi-stationary distance, with objects at larger distances requiring larger masses.

In conclusion, two possibilities for the source of CM Dra's timing variations remain valid: a mass of a few Jupiters on a two decade-long orbit, or an object on a century-to-millenium long orbit, with masses between 1.5 Jupiters and that of a very low mass star. Continued observations of the timing of CM Dra's eclipses over the next 5-10 years should, however, be decisive regarding the continued viability of the sinusoidal-fit model, and hence, about the validity of a Jovian-type planet in a circumbinary orbit around the CM Dra system.

Acknowledgements
Some of the observations published in this article were made with the IAC80 telescope operated by the Instituto de Astrofísica de Tenerife in the Observatorio del Teide, and with the INT telelescope operated by the Isaac Newtown Group of Telescopes in the Observatorio del Roque de los Muchachos. This research was supported by Grant ESP2004-03855-C03-03 of the Spanish Education and Science Ministry. Some material presented here is based on work supported by NASA under the grant NNG05GJ29G, issued through the Origins of Solar Systems Program.

References

 
Copyright ESO 2008