4 Statistical analysis

In the following section we perform a thorough statistical examination of the complete RV data set of the CES Long Camera survey. We test each star for variability, linear and curved RV slopes and significant periodic signals.

4.1 Probing for variability

To test if a star in the sample is variable we first apply the F-test ( $F = \sigma_{1}^{2} / \sigma_{2}^{2}$). We start by comparing the total variance of the RV results of one star ( $\sigma_{1}^{2}$) with the rest of the sample (P(F)₁ in Table 4). Since, for the small CES bandwidth, the RV precision and hence the scatter is a function of spectral type, we compare the stars with the mean scatter ( $\sigma_{2}^{2}$) of the remaining stars within the same spectral type bin (as described in Sect. 3).

For a second F-test, we take as the $\sigma_{2}^{2}$ the individual mean internal RV error and determine the P(F) (P(F)₂ in Table 4). With this we test if the star is more variable (or more constant) than its internal RV error would suggest.

Finally, we fit a constant to the RV results (with the zero-point as free parameter) and use the $\chi ^{2}$-statistic to check whether the constant model delivers a good description of the data. Again a small value of  $P_{\chi }(\chi ^{2})$means that the RV results for a star have a larger scatter than expected from their internal errors, which then can be interpreted as a sign for variability.

For the binaries ($\kappa$ For, HR 2400, HR 3677, $\alpha$ Cen A & B) the RV residuals after subtraction of either the known binary orbit ($\alpha$ Cen A & B) or the preliminary orbit we have found are taken for the variability tests. For $\alpha$ Cen A & B with the known binary orbit, HR 2400 with its linear trend and HR 3677 with its curved trend, $\nu$ the degrees of freedom equals N-1 (the velocity zero-point is adjusted), so that the F-test results for these stars are strictly valid only under the assumption that the binary orbit and the velocity trends are precisely known. In the cases of $\kappa$ For and $\iota$ Hor we take $\nu=N-6$. As criterion for a significant result we adopt the $99\%$-level (i.e. P < 0.01). Table 4 summarizes the results of these variability tests for all survey stars.

 

 

Table 4: F-test and $\chi ^{2}$-test results for the complete sample. P(F)₁ is the result of the F-test comparing the RV scatter with the average scatter within a spectral type bin, while P(F)₂ gives the results of the comparison with the mean internal RV error.  $\chi ^{2}_{\rm red}$ is the best (reduced) $\chi ^{2}$-value after fitting a constant to the RV data and $P_{\chi }(\chi ^{2})$ the corresponding probability.

Star	P(F)1	P(F)2	$\chi ^{2}_{\rm red}$	$P_{\chi }(\chi ^{2})$
$\zeta$ Tuc	0.023	0.073	1.92	0.0001
$\beta$ Hyi	0.13	0.12	1.51	3.5E-05
HR 209	0.5	0.35	0.98	0.50
$\nu$ Phe	3.6E-05	0.367	1.17	0.18
HR 448	0.34	0.40	0.70	0.85
HR 506	0.23	0.90	2.02	0.003
$\tau$ Cet	6.6E-11	0.016	0.73	0.98
$\kappa$ For $_{\rm Res}$	0.03	0.83	1.09	0.33
HR 753	0.63	0.203	0.32	0.90
$\iota$ Hor $_{\rm Res}$	0.009	4.9E-19	4.8	5.6E-45
$\alpha$ For	1.E-09	0.0014	1.82	6.9E-05
$\zeta ^{1}$Ret	0.57	0.71	1.55	0.09
$\zeta ^{2}$ Ret	0.67	0.05	1.69	0.0009
$\epsilon$ Eri	0.32	0.0059	1.86	3.5E-05
$\delta$ Eri	0.06	0.22	1.16	0.21
$\alpha$ Men	3.5E-06	0.37	0.90	0.64
HR 2400 $_{\rm Res}$	0.13	0.91	1.51	0.012
HR 2667	0.06	0.04	0.68	0.67
HR 3259	0.14	0.45	2.2	7.2E-05
HR 3677 $_{\rm Res}$	0.87	0.32	0.94	0.57
HR 4523	0.09	0.85	1.57	0.03
HR 4979	0.004	0.41	1.18	0.18
$\alpha$ Cen A $_{\rm Res}$	0.0001	0.78	0.70	0.94
$\alpha$ Cen B $_{\rm Res}$	0.88	0.13	1.26	0.12
HR 5568	0.001	0.002	0.38	0.99
HR 6416	0.09	0.0001	3.63	1.4E-18
HR 6998	0.70	0.28	0.55	0.99
HR 7373	0.02	0.83	0.77	0.62
HR 7703	0.65	0.75	0.81	0.76
$\phi ^{2}$ Pav	0.05	0.25	1.59	0.0003
HR 8323	0.85	0.56	1.02	0.43
$\epsilon$ Ind	0.39	0.0094	1.73	0.0001
HR 8501	3.1E-05	2.7E-07	2.26	3.7E-08
HR 8883	6.8E-10	0.005	4.21	9.3E-14
Barnard	0.0007	0.29	0.57	0.95
GJ 433	0.19	0.46	0.80	0.68
Prox Cen	2.2E-11	0.14	1.61	0.001

For 5 stars in the CES sample all 3 tests yielded P<0.01 and are thus clearly identified as RV variables: $\iota$ Hor (residuals), $\alpha$ For, HR 6416, HR 8501 and HR 8883. For $\iota$ Hor the cause of the residual variability was already discussed (stellar activity and possible second planet), while in the cases of $\alpha$ For ( $L_{\rm X}=525 \times 10^{27}~{\rm erg/s}$) and HR 8883 ( $L_{\rm X}=34175 \times 10^{27}~{\rm erg/s}$) also a high level of stellar activity appears to be responsible for the detected variability. We found strong Ca II H&K emission for HR 8883 by taking a spectrum with the FEROS instrument and the 1.5 m telescope on La Silla. We will show later that we can also determine the cause of variability for HR 6416 and HR 8501 (we will examine the presence of linear and curved trends in the RV data). The case of $\alpha$ For will be further discussed in more detail.

For 17 stars (HR 209, HR 448, $\kappa$ For, HR 753, $\zeta ^{1}$ Ret, $\delta$ Eri, HR 2400 (residuals), HR 2667, HR 3677 (residuals), HR 4523, HR 4979, $\alpha$ Cen B (residuals), HR 6998, HR 7373, HR 7703, HR 8323 and GJ 433) all 3 tests resulted in P>0.01 and thus no sign of variability whatsoever can be found for these stars. 5 stars are less variable than the rest of the sample (which of course also results in a low P(F)₁): $\nu$ Phe, $\tau$ Cet, $\alpha$ Cen A (residuals), HR 5568 and Barnard's star.

$\epsilon$ Eri and $\epsilon$ Ind both have a P(F)₁>0.01 (i.e. their overall RV scatter is quite typical for the CES sample) but a low P(F)₂ and $P_{\chi }(\chi ^{2})$, making them candidates for low-amplitude variations.

$\zeta$ Tuc, $\beta$ Hyi, HR 506, $\zeta ^{2}$ Ret and $\phi ^{2}$ Pav all appear unsuspicious in both F-tests but their RV data is not well fit by a constant function. Maybe an inclined linear or curved trend will describe these RV data better.

Finally, Prox Cen is more variable than the rest of the CES M-dwarf sample (of same spectral type) but not if compared to its own intrinsic mean RV error, however, for this star the  $P_{\chi }(\chi ^{2})$ for a constant fit is still less than 0.01.

4.2 RV trends: Linear and curved slopes

As the next step in the statistical analysis we examine whether the RV data of those stars which showed up as variable by the previous tests can be better described by a linear slope or a curved trend. For this purpose we determine the best-fit linear and parabolic function by $\chi ^{2}$-minimization and again compute  $P_{\chi }(\chi ^{2})$ for the following 12 stars: $\zeta$ Tuc, $\beta$ Hyi, HR 506, $\alpha$ For, $\zeta ^{2}$ Ret, $\epsilon$ Eri, HR 6416, $\phi ^{2}$ Pav, HR 8501, HR 8883, $\epsilon$ Ind and Prox Cen.

Table 5 shows the results of the slope and curvature tests. For 7 stars of the sample ($\zeta$ Tuc, HR 506, $\zeta ^{2}$ Ret, $\epsilon$ Eri, $\phi ^{2}$ Pav, HR 8883, Prox Cen) the resulting probabilities  $P_{\chi }(\chi ^{2})$ remain below 0.01. For these stars no linear or curved trend delivers a satisfactory explanation for their RV variations.

For $\epsilon$ Eri the $P_{\chi }(\chi ^{2})$-value for a curved trend is a magnitude higher than for a linear trend. As described in detail in Hatzes et al. (2000) there is evidence for a long-period ( $P\approx 2500$ days), low amplitude ( $K\approx19~{\rm m~s}^{-1}$) planet orbiting $\epsilon$ Eri. The difference we find between the linear and curved trend fits (see Fig. 12) might indicate the presence of the RV signature of this planet. The rms-scatter around the curved trend is with $12.9~{\rm m~s}^{-1}$ still larger than the mean internal error of $9.7~{\rm m~s}^{-1}$, which is not surprising for the modestly active star $\epsilon$ Eri.

Figure 12: Best-fit curved RV trend for $\epsilon$ Eri. With $\chi ^{2}_{\rm red}=1.7$ and $P_{\chi }(\chi ^{2})_{\rm curv}=0.0005$ this fit is still a poor description of the RV results. However, $P_{\chi }(\chi ^{2})_{\rm curv}$ is a magnitude higher than for a linear fit. The presence of the highly eccentric Keplerian signal of the planet described in Hatzes et al. (2000) is already indicated here.

For 5 stars ($\beta$ Hyi, $\alpha$ For, HR 6416, $\epsilon$ Ind and HR 8501) the linear slope tests resulted in a significant increase in  $P_{\chi }(\chi ^{2})$. This is a true indication for the presence of a linear trend in their RV behavior. However, for none of them, the probability of a curved trend is found to be significantly higher than for a linear one, except for $\beta$ Hyi where the difference is 5%. Table 6 gives the values of RV shift per day of the found linear trends and the remaining RV scatter around those trends. With the exception of $\alpha$ For the RV trends are all positive, with HR 8501 having the strongest slope of $+0.057~{\rm m~s}^{-1}~{\rm d}^{-1}$and $\epsilon$ Ind with the smallest variation of $+0.011~{\rm m~s}^{-1}~{\rm d}^{-1}$. Figures 13 to 17 display the best-fit trends in the RV results for these 5 stars.

  

Table 5: Results of linear slope and curvature tests. Probabilites $P_{\chi }(\chi ^{2})$ and $\chi ^{2}_{\rm red}$-values are given for the best-fit linear and parabolic functions. The first column ( $P_{\chi }(\chi ^{2})_{\rm constant}$) is repeated from Table 4 for comparison purposes.

Since HR 6416 and HR 8501 are known binary stars, the detected linear RV trend can be attributed to the binary orbital motion. The Hipparcos catalogue gives an angular separation of 8.658 arcsec for the HR 6416 binary. At a distance of 8.79 pc this implies a minimum separation of 76.1 AU. In the Gliese catalogue of nearby stars HR 6416 is listed as GJ 666A and the secondary GJ 666B as an M0V dwarf. We adopt $0.89~M_{\odot}$ as mass value for the G8V primary and for the secondary $0.52~M_{\odot}$ (after Gray 1988). Assuming a circular orbit and the minimum separation as the true separation and using Kepler's third law we can find an estimate for the RV acceleration for HR 6416. The orbital period is $\approx$565 yrs and the RV semi-amplitude $K \approx 1442~{\rm m~s}^{-1}$. Since the observing time span is not even a 1/100 of one orbital cycle we lineary interpolate to find an average acceleration of $\approx$ $0.028~{\rm m~s}^{-1}~{\rm d}^{-1}$ which is in good agreement with the detected trend of $0.032\pm0.003~{\rm m~s}^{-1}~{\rm d}^{-1}$.

The angular separation for the HR 8501 binary (GJ 853A & B) is given by the Gliese catalogue as 3.4 arcsec. This transforms into a minimum orbital separation of 46.3 AU at the distance of 13.61 pc. The spectral type of the secondary (GJ 853B) is unknown and therefore no mass estimate is possible. Assuming an M0V companion the average linear RV trend for the G3V primary ( $1.04~M_{\odot}$) would be $\approx$ $0.075~{\rm m~s}^{-1}~{\rm d}^{-1}$, slightly larger than the found trend of $+0.057\pm0.006~{\rm m~s}^{-1}~{\rm d}^{-1}$. The difference can easily be explained either by a larger orbital separation than the projected minimum value, the $\sin i$-effect, a lower mass of the secondary, the orbital phase corresponds to a steeper part of the RV-curve, or a combination of all these effects.

Figure 13: $\beta$ Hyi RV data and the best-fit linear slope (dashed line) and curved trend (dotted line). The rms scatter around these trends is $19.3~{\rm m~s}^{-1}$ and $19.5~{\rm m~s}^{-1}$ respectively. The found linear trend is  $+0.02\pm0.0024~{\rm m~s}^{-1}~{\rm d}^{-1}$.

Figure 14: Linear trend in the RVs of $\alpha$ For. The slope of $-0.04\pm0.007~{\rm m~s}^{-1}~{\rm d}^{-1}$ has a residual scatter of $50.9~~{\rm m~s}^{-1}$.

Figure 15: Best-fit linear trend for HR 6416, the RV shift per day is $+0.032\pm0.003~{\rm m~s}^{-1}$ and the rms scatter around this slope is  $19.4~{\rm m~s}^{-1}$.

Figure 16: Best-fit linear trend for $\epsilon$ Ind, the RV shift per day is $+0.012\pm0.002~{\rm ms}^{-1}$ and the rms scatter around this slope is $11.6~{\rm m~s}^{-1}$.

Figure 17: Best-fit linear trend for HR 8501, the RV shift per day is $+0.057\pm0.006~{\rm m~s}^{-1}$ and the rms scatter around this slope is  $23.4~{\rm m~s}^{-1}$.

 

 

Table 6: Linear RV trends and the rms scatter around these slopes. Uncertainties of the slopes correspond to the $\chi ^{2}=\chi ^{2}_{\rm best}+1$ range.

Star		RV trend	rms
		[ ${\rm m~s}^{-1}~{\rm d}^{-1}$]	[ ${\rm m~s}^{-1}$]
$\beta$ Hyi		$+0.02 ~~~ \pm 0.0024$	19.3
$\alpha$ For		$-0.04 ~~~ \pm 0.007$	50.9
HR 6416		$+0.032 \pm 0.003$	19.4
$\epsilon$ Ind		$+0.012 \pm 0.002$	11.6
HR 8501		$+0.057 \pm 0.005$	23.4

$\alpha$ For is also a known binary with an angular separation of 4.461 arcsec (from Hipparcos catalogue). At a distance of 14.11 pc this means $a\approx63$ AU. Adopting a mass for the F8V primary of $1.2~M_{\odot}$ (after Gray 1988) and assuming again an M0V secondary (the spectral type of the companion is unknown) we find an average acceleration of  $0.04~{\rm m~s}^{-1}~{\rm d}^{-1}$ which is exactly the value of the RV trend we observe for $\alpha$ For. However, the rms scatter around this trend is large with $\approx$ $51~{\rm m~s}^{-1}$. Most of this residual scatter is probably caused by the high stellar activity of $\alpha$ For. Another possible explanation for an increase of the RV scatter can be contamination of the spectra of the primary by the secondary. Due to field rotation at the Nasmyth focus of the CAT telescope, light from the close secondary also entered the CES during some observations when the slit was aligned with the binary axis. This could at least to some degree have contributed to the larger scatter (the same is true for HR 8501 with an even smaller angular separation).

Neither $\beta$ Hyi nor $\epsilon$ Ind are known binary stars. The detected linear RV trends are thus caused by previously unknown companions. The linearity of both trends points towards distant stellar companions. However, both stars also represent candidates for having very long-period (P>20 yrs) planetary companions. Follow up observations using the upgraded CES, now equipped with the Very Long Camera (VLC), and the 3.6 m telescope were already performed and are still in progress. Analysis of the new data and the combination with the RV results presented here will show whether the linearity of both trends continues.

4.3 Period search

To search for periodic signals in the complete RV results of the CES sample we again use the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). We estimate the False Alarm Probability (FAP) of a peak in the power spectrum by employing a bootstrap randomization method. In this bootstrap approach the actual RV measurements are randomly redistributed while keeping the times of observations fixed (Kürster et al. 1997; Murdoch et al. 1993). The major advantage of this method is the fact that the FAP-levels can be derived without any assumptions on the underlying noise distribution (like e.g. a Gaussian). We will continue to follow this bootstrap philosophy also in determining upper mass-limits for planets in Sect. 5.

The search interval is 2 to 5000 days. For each star we perform 10 000 bootstrap randomization runs to estimate the FAP of the maximum peak in the power spectrum. Since the validity of the FAP resulting from random redistribution is lower in the case of high temporal concentration at one or several points (i.e. "data clumping'' when during one night a large number of measurements were taken while in other nights this number is significantly lower) we also perform the analysis on the RV set binned in nightly averages.

Table 7 summarizes the periodogram results for all survey stars (in the cases of $\kappa$ For, HR 2400, HR 3677 and $\alpha$ Cen A & B the period search was performed on the RV residuals after subtraction of the binary orbit).

  

Table 7: Periodogram results. Power₁, Period₁ and FAP₁ are the results of the original RV data set (N₁), while Power₂, Period₂ and FAP₂ gives the values for the nightly averages (N₂).

The periodogram analysis of the original RV data found periods with FAP₁ (for the original un-binned RV data set) below 0.001 at: $\beta$ Hyi, $\tau$ Cet, $\iota$ Hor, $\alpha$ Cen A & B residuals, $\phi ^{2}$ Pav and HR 8501. However, none of them, with the exception of $\iota$ Hor and HR 8501, reveal a significant signal in their nightly averaged RV data set. In the cases of $\beta$ Hyi and $\alpha$ Cen A & B the significance of the low FAP₁ is decreased by high temporal concentration of their RV data (the periodogram results for $\alpha$ Cen A & B are discussed in detail in Endl et al. 2001a). The significant P=5000 day signal for HR 8501 is caused by the found linear RV trend (see Fig. 17), since 5000 days is the maximum period searched for. The same happens in the case of HR 6416, although the P=5000 day signal is no longer significantly recovered in the nightly averaged data set.

Thus $\iota$ Hor remains the single case in the CES Long Camera survey which shows a convincing periodic RV variation of statistical significance due to an orbiting planet.

The results for $\phi ^{2}$ Pav are interesting, since the FAP₂ of the P=7 day signal is still low with 0.003 and just marginally above the 1 permill threshold. We will discuss the case of $\phi ^{2}$ Pav later.

4.4 Activity induced RV excess scatter

Several studies (Saar & Donahue 1997; Saar et al. 1998; Santos et al. 2000; Paulson et al. 2002) have investigated and discussed the relationship between stellar activity and RV scatter induced by activity related phenomena like surface inhomogenities (spots) and variable granulation pattern. Saar et al. (1998) found a correlation between an excess RV scatter (exceeding the expected scatter due to internal measurement uncertainties) and the rotational speed of the G and K dwarfs in the sample of the Lick planet search program. Interestingly, Paulson et al. (2002) measured simultaneously the Ca II H&K emission and RVs of members of the Hyades cluster and showed that only for a few stars (5 out of 82) the chromospheric activity is correlated with an RV excess.

In the case of the CES Long Camera survey we can estimate $P_{\rm Rot}$ for 13 stars, which reveal an excess RV scatter compared to their internal errors, and examine whether this excess scatter is correlated with intrinsic stellar activity. Stellar rotational periods are either taken from Saar & Osten (1997) or estimated using Eqs. (3) and (4) in Noyes et al. (1984) and the Ca II H&K results from Henry et al. (1996). Figure 18 displays the RV excess scatter plotted vs. $P_{\rm Rot}$, which shows a general increase in the excess scatter with shorter rotational periods. The sample of stars included in the diagram consist of 2 F-type stars ($\zeta$ Tuc, $\nu$ Phe), 7 G-type stars ($\iota$ Hor, $\zeta ^{2}$ Ret, HR 4523, HR 4979, HR 6416, HR 8323, HR 8501) and 4 K-type stars ($\epsilon$ Eri, $\delta$ Eri, $\alpha$ Cen B, $\epsilon$ Ind). Although we do not have $P_{\rm Rot}$-values for $\alpha$ For and HR 8883, the observed large excess scatter for these two stars is probably also due to fast rotation and enhanced stellar activity. As already mentioned, both stars are bright X-ray sources and HR 8883 displays prominent Ca II H&K emission. This result confirms the general picture of increased RV-jitter due to enhanced stellar activity.

Figure 18: Correlation between RV excess scatter and stellar rotational period. F-type stars are shown as circles, G-type stars as full diamonds and K-type as triangles. The excess scatter increases with decreasing  $P_{\rm Rot}$ (with $\iota$ Hor having the highest value of excess scatter).