# The Lick Planet Search : Detectability and Mass Thresholds

## 1 Introduction

In the past few years, high precision radial velocity surveys have had remarkable success in the discovery of planetary-mass companions around nearby solar-type stars (for reviews, see Marcy & Butler 1998 and Marcy, Cochran & Mayor 1999). Searches for companions (Campbell et al. 1988; McMillan et al. 1994; Mayor & Queloz 1995; Walker et al. 1995; Cochran et al. 1997; Noyes et al. 1997; Marcy & Butler 1992, 1998) have been carried out with Doppler velocity precision $`10\mathrm{m}\mathrm{s}^1`$, although $`3\mathrm{m}\mathrm{s}^1`$ has been achieved at Lick observatory for chromospherically quiet stars (Butler et al. 1996). There are now 17 companions known with masses (the observable is $`M_p\mathrm{sin}i`$, where $`i`$ is the angle of inclination of the orbit with respect to the line of sight) below 10 Jupiter masses. In total, several hundred stars have been monitored for timescales of 3 years to more than 11 years. The detections so far suggest that a few percent of solar type stars harbor companions of a Jupiter mass or more within a few AU.

These objects have raised many questions regarding the distribution of the mass and orbital radius of planetary-mass companions, and the relation of these systems to our own solar system and its giant planets. For example, a surprise was the discovery of Jupiter-mass companions in close proximity to their host star. Of the 17 companions within $`2.5\mathrm{AU}`$, 13 have semi-major axis $`a<0.5\mathrm{AU}`$ and five have $`a<0.1\mathrm{AU}`$. The archetypal example is the companion orbiting 51 Pegasi (Mayor & Queloz 1995) which has a mass ($`M_p\mathrm{sin}i`$) of 0.44 Jupiter masses ($`M_J`$) and an orbital radius $`a=0.05\mathrm{AU}`$, eight times closer than Mercury’s orbit about the Sun. The orbital parameters of the 17 planetary-mass companions are listed in Table 1.

Unfortunately, gleaning the true distribution of companions is complicated by selection effects which favor the detection of massive, close companions. It is necessary to establish detection thresholds for searches for planetary-mass companions before the observations can be fully interpreted. Walker et al. (1995) monitored 21 bright solar type stars for 12 years. They carried out a detailed statistical analysis, and from their upper limits could exclude companions with $`M_p\mathrm{sin}i1`$$`3M_J`$ for periods less than the duration of their observations ($`12`$ years). Nelson & Angel (1998) used a simple analytic formalism, together with comparisons with real data, to investigate the dependence of detection thresholds on the number and duration of observations and the Doppler errors.

Our aim in this paper is to place the confirmed companions from the Lick radial velocity survey in context by an analysis of the null detections. The ongoing Lick survey consists of more than eleven years of precision Doppler velocity measurements of 107 stars (a few of which have been added or dropped as the survey progressed), and 200 new stars have recently been added (Fischer et al. 1999). The original program has so far identified six planetary-mass companions (70 Vir, 47 UMa, 55$`\rho `$ Cnc, $`\tau `$ Boo, $`v`$ And, GL876, they are marked in Table 1), codiscovered the companion to 16 Cyg B, and confirmed two discoveries made by other groups (51 Peg and $`\rho `$ CrB, see Table 1 for orbital parameters and references).

We search for periodicities and place upper limits using the “floating-mean periodogram”, an extension of the well-known Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) in which we fit sinusoids to the data, but allow the zero-point of the sinusoid to “float” during the fit. This approach was adopted by Walker et al. (1995), who were interested in obtaining correct upper limits for periods greater than the duration of the observations. We show here that allowing the mean to float is crucial to account for statistical fluctuations in the mean of a sampled sinusoid. The traditional Lomb-Scargle periodogram fails in precisely the regime where we demand it be robust, namely when the number of observations is small, the sampling is uneven or the period of the sinusoid is comparable to or greater than the duration of the observations. We carefully consider the correct normalization of the periodogram, an issue which has been of some debate in the literature (Horne & Baliunas 1986; Gilliland & Baliunas 1987; Schwarzenberg-Czerny 1996).

The plan of the paper is as follows. We begin in §2 by describing the observations and estimating the velocity variability we expect to see due to measurement error and intrinsic stellar effects. We discuss in §3 our methods for searching for the signatures of companions in the radial velocity data. We look for variability in excess of our prediction, long term trends and periodicities. We list those stars whose data show interesting variability or periodicities which may indicate the presence of a yet unconfirmed companion. In §4, we place upper limits on the mass of a possible companion as a function of orbital radius for each star in the sample. We continue in §5 with a discussion of the implications of our results for the occurrence rate of planetary-mass companions to solar type stars, and their distribution in mass and orbital radius. We present our conclusions in §6. In Appendix A, we give a brief derivation of the periodogram, and in Appendix B, we discuss its normalization.

## 2 The Observations

The Lick radial velocity survey is now more than eleven years old (Marcy & Butler 1992, 1998). Precise radial velocity measurements (current precision $`5\mathrm{m}\mathrm{s}^1`$, Butler et al. 1996) are made with the Lick 3m telescope by using an echelle spectrograph and a comparison iodine reference spectrum. The exposure time is ten minutes for a star with $`V=5`$, allowing several observations per year for each star in the sample. In this paper, we present an analysis of observations of 76 F, G, and K type stars in the original survey<sup>1</sup><sup>1</sup>1There have been observations of 29 M dwarfs as part of the survey, but we do not include these in our analysis as they are faint ($`V>7`$) and suffer from large measurement uncertainties. These stars are part of a sample being monitored with the Keck telescope. An analysis of their radial velocity variability will be presented elsewhere.. A summary of the observations is given in Table 2. For each star, we give its HR and HD catalog number, spectral type and rotation period. We list the number and duration of the observations, the typical internal velocity error and the rms scatter of the data. Radial velocities are available upon request from G. Marcy.

### 2.1 The Distribution of Errors and Intrinsic Stellar Variability

Two sources of variability can mask velocity variations due to a companion: measurement uncertainties and intrinsic stellar variability. In this section, we attempt to quantify these effects.

The uncertainty in the radial velocity measurement $`v`$ is estimated for each observation from the dispersion of the velocities measured by different spectral segments of the spectrometer. An upgrade to the spectrograph optics and improvements in modeling in November 1994 led to an increase in the Doppler precision from $`\sigma _D10`$$`15\mathrm{m}\mathrm{s}^1`$ to $`\sigma _D5\mathrm{m}\mathrm{s}^1`$, with $`\sigma _D3\mathrm{m}\mathrm{s}^1`$ achievable in the best cases. This improved Doppler precision is dominated by photon statistics (for a detailed discussion, see Butler et al. 1996). In this paper, we shall refer to the pre-November 1994 data as “pre-fix” and post-November 1994 data as “post-fix”. Table 2 gives the average internal error $`<\sigma _D>`$ before and after November 1994 for each star.

Intrinsic stellar variability arises from magnetic activity or rotation of features across the stellar surface, such as sunspots or inhomogeneous convective patterns (Saar & Donahue 1997). Saar, Butler & Marcy (1998, hereafter SBM98) used the Lick radial velocity variations (post-fix data only) to characterize the relationship between the rotation period, $`P_{\mathrm{rot}}`$, of a star and its intrinsic velocity variability, $`\sigma _V`$. They found that the variability in excess of the internal errors could be explained by simple models of sunspot rotation and inhomogeneous convective flows. We use their results to estimate the typical intrinsic variability associated with such effects as a function of stellar rotation period<sup>2</sup><sup>2</sup>2 At first sight, it may seem circular to use the work of SBM98 which was based on the Lick data set (the post-fix data only). However, our approach is to use their results to characterize the average variability typical for a star in the survey with a particular rotation period.. After inspecting Figure 2 of SBM98, we find $`\sigma _V=10\mathrm{m}\mathrm{s}^1(12\mathrm{days}/P_{\mathrm{rot}})^{1.1}`$ for G and K type stars and $`\sigma _V=10\mathrm{m}\mathrm{s}^1(10\mathrm{days}/P_{\mathrm{rot}})^{1.3}`$ for F stars. The rotation period for each star (taken from Soderblom 1985, Baliunas, Sokoloff & Soon 1996, or Fischer 1999, private communication) is given in Table 2.

We obtain the total estimated variability for each data point by adding the intrinsic variability to the internal error in quadrature, $`\sigma _j^2=\sigma _V^2+\sigma _{D,j}^2`$ (here we label each data point with the index $`j`$). How well does this estimate reproduce the scatter in the data? After a preliminary analysis of the post-fix data (using the methods of §§3.1 and 3.2), we selected a subset of 26 stars that had no excess variability or long term trends. In Figure 1, we plot a histogram (solid line) of individual radial velocity measurements $`v_j`$ divided by the estimated variability $`\sigma _j`$ for this subset of stars. For each star, we have subtracted the weighted mean of the velocities. The upper panel shows the pre-fix data; the lower panel shows the post-fix data. The dotted histogram in each case shows a Gaussian distribution with unit variance. If the scatter in the velocities were Gaussian with variance $`\sigma _j^2`$ for each point, the dotted and solid histograms would match. They do not, indicating more scatter in the data than we expect. In addition, the pre-fix and post-fix $`v/\sigma `$ distributions are different.

The excess scatter may be due to a combination of underestimated internal errors and systematic errors in the velocities, particularly for the early observations. We have chosen to augment the internal errors by multiplying by a constant factor to force the observed $`v/\sigma `$ distribution to match a Gaussian with unit variance. In this way, we bring the pre-fix and post-fix $`v/\sigma `$ distributions into agreement, and we are confident that we have not underestimated the errors. For all 76 stars in our sample, we multiply the pre-fix internal errors by 1.7 and the post-fix internal errors by 1.4. The dashed histograms in Figure 1 show the $`v/\sigma `$ distributions using these rescaled internal errors. The distribution is unchanged if we remove the 10 chromospherically-active stars in this subset which have $`P_{\mathrm{rot}}12\mathrm{days}`$. Hereafter, we refer to the augmented internal errors as simply “internal errors”.

In Figure 2, we show the effect of including the intrinsic variability prediction $`\sigma _V`$ for stars of different rotation period. For each star in the subset of 26 stars of Figure 1, we plot $`\chi _\nu ^2=(v_j/\sigma _j)^2/(N1)`$ as a function of $`P_{\mathrm{rot}}`$. We show $`\chi _\nu ^2`$ evaluated using the internal errors only (crosses) and including the intrinsic variability added in quadrature (squares). The extra variability of stars with short rotation periods ($`P_{\mathrm{rot}}12`$ days) is shown by the large uncorrected $`\chi ^2`$ values for these stars. The mean value of $`\chi _\nu ^2`$ in Figure 2 is less than unity, indicating that our procedure may have overestimated the internal errors somewhat. However, we prefer to err on the side of overestimation.

We use the variability estimate $`\sigma _j`$ in two ways. The first is to identify those stars which show much more variability than we might expect given their rotation periods (§3.2). The second is to weight the data points when we look for periodicities (§3.4). The large difference between the pre-fix and post-fix data makes it important to give less weight to the early data points. Not only are the pre-fix errors larger, they are less well characterized than the post-fix errors. One might question the value of including the low quality early data points at all. However, they are important because they extend the time baseline, allowing us to search for longer period companions.

## 3 Search for Companions

The velocity amplitude $`K`$ of a star of mass $`M_{}`$ due to a companion with mass $`M_p\mathrm{sin}i`$ with orbital period $`P`$ and eccentricity $`e`$ is

$$K=\left(\frac{2\pi G}{P}\right)^{1/3}\frac{M_p\mathrm{sin}i}{(M_p+M_{})^{2/3}}\frac{1}{\sqrt{1e^2}}.$$
(1)   
For a circular orbit with $`M_pM_{}`$, the velocity variations are sinusoidal with amplitude

$$K=28.4\mathrm{m}\mathrm{s}^1\left(\frac{1\mathrm{year}}{P}\right)^{1/3}\left(\frac{M_p\mathrm{sin}i}{M_J}\right)\left(\frac{M_{}}{M_{}}\right)^{2/3},$$
(2)   
where $`M_J`$ is the mass of Jupiter. The orbital period is related to the orbital radius by Kepler’s law,

$$P=1\mathrm{year}(a/\mathrm{AU})^{3/2}(M_{}/M_{})^{1/2}.$$
(3)   
For example, the companion to 51 Peg ($`a=0.05`$ AU, $`M_p\mathrm{sin}i=0.44`$) induces a velocity amplitude $`K=56\mathrm{m}\mathrm{s}^1`$, whereas Jupiter ($`a=`$5.2 AU, $`P=11.9`$ years) gives $`K=12.5\mathrm{m}\mathrm{s}^1`$ for the Sun.

In this section, we describe the methods we use to search for such a signal, and present our results. For each star in Table 2, we first ask if there is a significant long term trend in the data, and if so we subtract it (§3.1). We then ask if the observed scatter in the data is consistent with the expected variability (§3.2). To search for periodicities (§3.4), we fit sinusoids to the data, employing a generalization of the well-known Lomb-Scargle periodogram (§3.3). By using sinusoids as our basic model for the data, we are strictly assuming circular orbits, although we find that the periodogram gives a good estimate of the orbital period even for eccentric orbits. This is because, to lowest order in the eccentricity, the radial velocity signal from an eccentric orbit has its main component at the orbital frequency (with smaller components at multiples of the orbital frequency). We discuss a possible extension of the periodogram to Kepler orbits in the conclusions (§6). The normalization of the periodogram has been of some question in the literature, so we discuss our choice of normalization in Appendix B. We close this section by summarizing our results and indicating those stars that show interesting variability or evidence for companions (§3.5).

### 3.1 Long Term Trends

We first ask if there is a significant long term trend, on a timescale much greater than the duration of the observations (i.e. $``$ 10 years). For each star, we fit a straight line $`v_j=at_j+b`$ to the measured velocities. When calculating $`\chi ^2`$ for the fit, we weight each point by the inverse square of the estimated error, $`w_j=1/\sigma _j^2`$. We give the best-fit slope and its uncertainty (as derived from the least-squares fit) in Table 3.

To assess the significance of each slope, we ask if the coefficient of the linear term is significantly non-zero. We use the F-test to compare the weighted sum of squares of residuals from the straight line fit $`\chi _{N2}^2`$ (two free parameters) to the weighted sum of squares about the mean $`\chi _{N1}^2`$ (one free parameter). If there is no long term trend and the residuals are Gaussian distributed, the statistic

$$F=(N2)\frac{\chi _{N1}^2\chi _{N2}^2}{\chi _{N2}^2}$$
(4)   
follows Fisher’s $`F`$ distribution (Bevington & Robinson 1992) with $`1`$ and $`N2`$ degrees of freedom. Given the distribution $`F_{1,N2}`$, we calculate the probability that $`F`$ would exceed the observed value $`F_{\mathrm{obs}}`$ purely by chance (the false alarm probability). We give the F-test false alarm probabilities for each star in Table 3.

The slopes listed in Table 3 contain much information about possible companions at long periods. For the purposes of this paper, however, we are interested in identifying slopes which would directly affect our search for companions with $`P30`$ years. Thus, we mark with $``$ in Table 3 those stars which have a slope $`5\mathrm{m}\mathrm{s}^1\mathrm{yr}^1`$ and an F-test false alarm probability $`<10^5`$. We adopt a higher threshold ($`10^5`$) than elsewhere in this paper because systematic effects in the pre-fix data or variations due to magnetic activity on long timescales can imitate a slope. The stars marked with $``$ in Table 3 are HR 219a, HR 753, HR 1325, HR 2047, HR 5544a, HR 5544b, HR 6623, HR 7672 and HR 8086<sup>3</sup><sup>3</sup>3HR 8085, the companion to HR 8086, would have exhibited a slope had we removed the effect of secular acceleration, see §3.5..

### 3.2 Excess Variability

We now ask if the scatter in the velocities is consistent with our predicted variability for each star (§2.1). The $`\chi ^2`$ about the mean, $`\chi _{N1}^2`$, or (for those stars marked $``$ in Table 3) about the best-fit straight line, $`\chi _{N2}^2`$, gives us a measure of variability. We use the $`\chi ^2`$ distribution with $`Nm`$ degrees of freedom to test if the velocities are consistent with being drawn from a Gaussian distribution with variance $`\sigma _j^2`$. Here $`m`$ is the number of parameters in the model of the data, $`m=2`$ for a straight line or $`m=1`$ for a mean. A small false alarm probability indicates there is more variability in the data than we expect. The results of this test are shown in Table 4. We indicate with a $``$ those stars which have false alarm probabilities less than 1%. We choose a 99% threshold so that there will be no more than one false signal in our sample of 76 stars. The stars which show significant variability are HR 88, HR 166, HR 2047, HR 4345, HR 5273, HR 5544b, HR 5553, HR 7061, GL 641, GL 688, and the stars with confirmed planetary-mass companions (listed in Table 1).

### 3.3 The Floating-Mean Periodogram

In this section, we test for periodicity in the data using what we refer to as the “floating-mean periodogram”, a generalization of the well-known Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). We define the floating-mean periodogram as follows. For each trial frequency $`\omega =2\pi /P`$, we start with a simple model for the data, namely a sinusoid plus a constant term,

$$f_j=A\mathrm{cos}\omega t_j+B\mathrm{sin}\omega t_j+C,$$
(5)   
where $`t_j`$ are the observation times. We perform a linear least squares fit of this model to the data, to determine the constants $`A`$,$`B`$ and $`C`$. The periodogram is an “inverted” plot of the $`\chi ^2`$ of this fit as a function of frequency. We define the floating-mean periodogram power $`z(\omega )`$ as

$$z(\omega )=\frac{(N3)}{2}\frac{\chi _{N1}^2\chi ^2(\omega )}{\chi ^2(\omega _0)},$$
(6)   
where $`\chi ^2(\omega )=w_j[v_jf_j(\omega )]^2`$ is the $`\chi ^2`$ of the fit, $`\omega _0`$ is the best-fit frequency (i.e. the frequency that gives the maximum periodogram power) and $`\chi _{N1}^2`$ is the weighted sum of squares about the mean. When calculating $`\chi ^2`$, we weight each data point by $`w_j1/\sigma _j^2`$ as in §3.1. We use a linear least squares fitting algorithm from Press et al. (1992) to fit equation (5) to the data and find $`\chi ^2(\omega )`$. The choice of normalization of the periodogram has been a subject of some debate in the literature; we discuss this in detail in Appendix B. We normalize by the weighted sum of squares of the residuals to the best fit sinusoid, $`\chi ^2(\omega _0)`$.

The Lomb-Scargle periodogram is obtained by considering a fit of a sinusoid only, i.e. the case $`C0`$ in equation (5) (we sketch the derivation of the Lomb-Scargle formula from the least-squares approach in Appendix A). This means that the zero-point of the sinusoid is assumed to be known already. In practice, the mean of the data is taken as an estimate of the zero-point and is subtracted from the data before applying the Lomb-Scargle formula. In our approach, the zero-point of the sinusoid is an additional free parameter at each frequency, i.e. the mean of the data is allowed to “float” during the fit.

This approach has been adopted by other authors. Ferraz-Mello (1981) was the first to do so, defining the “date-compensated discrete Fourier transform”. Walker et al. (1995) generalized to the case where a straight line or quadratic function was subtracted from the data, defining a “correlated periodogram”. Most recently, Nelson & Angel (1998) included a constant term in their Monte Carlo experiments. These authors were concerned about the suppression of periodogram power for periods greater than the duration of the observations. We show here that allowing the mean to float is important under much more general circumstances.

We now provide some examples which show that allowing the mean to float is crucial if the number of observations is small, the sampling is uneven or there is a period comparable to the duration of the observations or longer. Figure 3 shows simulated data of a companion in a circular orbit with $`P=9.6`$ years and $`K=15\mathrm{m}\mathrm{s}^1`$ ($`a=4.5\mathrm{AU}`$, $`M_p\mathrm{sin}i=1.1M_J`$). We use the observation times and errors for one of the stars in our sample, HR 222. The upper panel shows the velocity measurements as a function of time. The dashed line shows the mean of the data, which is about $`10\mathrm{m}\mathrm{s}^1`$ greater than the correct zero-point of the sine wave. The solid line shows the best-fit sinusoid when the mean is allowed to float. The traditional and floating-mean periodograms are shown in the lower panel. The power at long periods is significantly less in the traditional periodogram than the floating mean periodogram. Black & Scargle (1982) first noted this effect in their analysis of astrometric data, where they showed that proper motions could significantly reduce the measured amplitude of long period signals (see Figure 3 of Black & Scargle 1982).

The upper panel in Figure 4 shows 20 velocities obtained at Lick for the star HR 5968, which has a known companion with $`K=67\mathrm{m}\mathrm{s}^1`$ and $`P=39.6`$ days (Noyes et al. 1997, Table 1). The duration of these observations is 1.2 years. By chance, most of the measurements lie above the zero-point of the orbit, giving a $`20\mathrm{m}\mathrm{s}^1`$ difference between the mean of the data and the correct zero-point. We plot the best-fit sinusoid with a fixed mean as a dashed line and with a floating mean as a solid line. The traditional periodogram does not detect a significant period; it identifies a period of $`43.9`$ days, but with false alarm probability 8%. The floating-mean periodogram gives a very significant detection at the correct period of $`40.0`$ days (see Table 5). Thus allowing the mean to float is not only important at long periods, but is crucial to account for statistical fluctuations when the number of observations is small. A similar situation could occur due to uneven sampling. If a sinusoid is sampled at nearly the same phase each cycle (e.g. for $`P1`$ year), the mean of the data could be significantly different from the zero-point of the sinusoid, and the periodicity thus go undetected.

If the number of data points is large and the periodicity is well-sampled, the mean of the data does give a good estimate of the correct zero-point. However, the traditional periodogram fails in precisely the regime where we require it to be robust; namely, when the number of observations is small, the duration of the observations is limited, or the sampling uneven. For this reason, we adopt the floating-mean periodogram, despite it being less computationally-efficient than the simple Lomb-Scargle formula.

Following Walker et al. (1995), we make a further generalization in the case where we subtract a straight line from the data (those stars marked $``$ in Table 3). In this case we must fit a sinusoid plus a straight line to the data at each frequency (i.e. add a term $`t_j`$ to eq. ). For these stars, we use $`\chi _{N2}^2`$ in place of $`\chi _{N1}^2`$ in our definition of the periodogram, and replace $`\chi ^2(\omega )`$ by the $`\chi ^2`$ from the straight line plus sinusoid fit. This gives a general formula for the floating-mean periodogram power (as Walker et al. 1995, eq. \[A2\])

$$z(\omega )=\frac{(Nm2)}{2}\frac{\chi _{Nm}^2\chi _{Nm2}^2(\omega )}{\chi _{Nm2}^2(\omega _0)},$$
(7)   
where $`m=1`$ for a mean or $`m=2`$ for a straight-line. Comparing equation (7) with equation (4), we see that the periodogram is similar to the F-statistic, measuring how much the fit is improved by introducing the two extra sinusoid parameters.

### 3.4 Application of the Periodogram

#### 3.4.1 Search for Periodicities

For each star, we plot the periodogram and look for the frequency which gives the maximum power $`z_{\mathrm{max}}`$. The periodogram power $`z`$ is a continuous function of frequency $`f=\omega /2\pi `$. However, the finite duration of the observations $`T`$ gives each periodogram peak a finite width $`\mathrm{\Delta }f1/T`$. Thus in a frequency range $`\mathrm{\Delta }f`$, there are roughly $`T\mathrm{\Delta }f`$ peaks. To make sure we sample all the peaks, we evaluate $`4T\mathrm{\Delta }f`$ periods between 2 days and 30 years<sup>4</sup><sup>4</sup>4The “average” Nyquist period of our observations ($`P_{\mathrm{Nyq}}2T/N`$) is a few months. However, the uneven sampling gives information on much shorter periods (perhaps much shorter than the minimum spacing between observations, see Eyer & Bartholdi 1998). The minimum period we investigate is 2 days. The maximum period of 30 years is a few times greater than the typical duration of the observations.. Monte Carlo tests indicate that this gives adequate sampling of the periodogram. We evaluate $`z`$ at evenly-spaced frequencies using a linear least squares fitting algorithm from Press et al. (1992).

To assess the significance of a possible detection, we test the null hypothesis that the data are pure noise. We ask, what is the false alarm probability associated with $`z_{\mathrm{max}}`$, or how often would noise fluctuations conspire to give a maximum power larger than that observed? We use Monte Carlo tests to find the false alarm probability<sup>5</sup><sup>5</sup>5The false alarm probability increases as the frequency range searched $`\mathrm{\Delta }f`$ increases (Schwarzenberg-Czerny 1996, 1997a, 1998 refers to this as the bandwidth penalty). As we discuss in Appendix B, although the distribution of $`z`$ at one particular frequency is easy to write down analytically, the number of independent frequencies in a frequency range $`\mathrm{\Delta }f`$ is not. For this reason, we adopt a Monte Carlo approach. This also allows us to check for non-Gaussian effects.. For each star, we make fake data sets of either a mean or straight line plus noise. We then perform the same analysis as for the original data. We find the maximum periodogram power, and ask, in what fraction of trials does $`z_{\mathrm{max}}`$ exceeds the observed value? This fraction is the false alarm probability.

We add noise to the simulated data sets in two ways. One is to add Gaussian deviates with the same variance as the observed velocities, keeping the same relative weights and observation times. The second is to take the observed velocities and randomize them, keeping the sample times fixed (the so-called bootstrap method, Press et al. 1992). In this case, we randomize the pre-fix and post-fix velocities separately, to account for their different $`v/\sigma `$ distributions. We find that both approaches give false alarm probabilities which are similar for most stars. We have also applied the analytic distribution given in Appendix B to our results (see Table 9), and fit for the number of independent frequencies $`M`$. The resulting analytic false alarm probabilities agree well with our Monte Carlo calculations.

The results are given in Table 5. For each star, we give the maximum power $`z_{\mathrm{max}}`$, the corresponding best-fit period and velocity amplitude, and the false alarm probabilities determined from 400 Monte Carlo trials for each star. We mark with $``$ those stars which show false alarm probabilities less than 1%. Apart from those stars with confirmed companions (Table 1), these are HR 509, HR 937, HR 996, HR 1084, HR 1729, HR 2047, HR 3951, HR 4345, HR 4496, HR 4540, HR 4983, HR 5019, HR 5273, HR 5544a, HR 5568, HR 6623, HR 7061, HR 7602, HR 7672, and HR 8086. Again, our motivation for choosing a 99% detection threshold is that we then expect no more than one false detection in our sample of 76 stars.

We find that 20 stars have more than one peak in the periodogram with false alarm probability $`<1\%`$. The search for multiple companions is beyond the scope of this paper, however, we carried out a simple test of whether these secondary peaks were due to aliasing of the primary period by the finite sampling. We subtracted off the best fit sinusoid from the data and looked at a periodogram of the residuals to see if the secondary peaks remained. Only in two cases did they remain, HR 458 which has a second peak at 1210 days, and HR 509 which has a second peak at 60 days.

#### 3.4.2 Variability

The average power $`\overline{z}`$ (averaged over frequency) is an additional indicator of variability. If the data are drawn from a Gaussian distribution, we expect the mean power to be $`\overline{z}1`$ (the mean value of the $`F`$-distribution). Again using a Monte Carlo method, we calculate the false alarm probability associated with the observed $`\overline{z}`$. Namely, we ask in what fraction of simulated data sets is the mean power larger than the observed value? A low false alarm probability in this test indicates either a periodicity may be present (the uneven sampling results in “leakage” which contaminates the background level), some non-Gaussian behavior or some kind of “broad-band” variability (for example as might be expected from magnetic activity). The results are shown in Table 6. Several stars show significantly high $`\overline{z}`$ at the 99% level. Apart from those stars with confirmed companions (Table 1), they are HR 88, HR 166, HR 509, HR 937, HR 996, HR 1084, HR 1729, HR 2047, HR 4345, HR 4496, HR 4540, HR 4983, HR 5019, HR 5273, HR 5544a, HR 5568, HR 5868, HR 5914, HR 5933, HR 6623, HR 7061, and HR 7602.

### 3.5 Results of our Search for Companions

Our sample of 76 stars contains eight stars with companions of planetary-mass, HR 3522, HR 4277, HR 458, HR 5072, HR 5185, HR 5968, HR 7504, and HR 8729. All of these companions give highly significant periodogram peaks. The orbital period and velocity amplitude obtained from the periodogram agree well with a Keplerian fit, although not surprisingly, the velocity amplitude is less than that of the Kepler fit for eccentric orbits. The orbital parameters for these planetary-mass companions and references are given in Table 1.

Five stars have companions with $`M_p\mathrm{sin}i>15M_J`$ which have been reported by previous authors. They are HR 2047 ($`M\mathrm{sin}i=0.15M_{}`$, $`P=14.2\mathrm{yr}`$; Irwin, Yang, & Walker 1992), HR 5273 ($`M\mathrm{sin}i=0.4M_{}`$, $`a=5\mathrm{AU}`$; Kamper 1987), HR 5553 ($`K20\mathrm{km}\mathrm{s}^1`$, $`P=125`$ d; Beavers & Salzer 1983), HR 6623 (Cochran & Hatzes 1994 report a slope of $`40\pm 5\mathrm{m}\mathrm{s}^1\mathrm{yr}^1`$) and GL688 ($`K=5.7\mathrm{km}\mathrm{s}^1`$, $`P=83.7`$ d; Tokovinin 1991). These massive companions explain the variability seen in these stars, and the trend in those two, HR 2047 and HR 6623, for which we have not yet seen a complete orbital period. They provide interesting examples of the limitations of the periodogram. There are 20 observations of HR 5553. The orbital period is correctly identified by the periodogram as 125 days, but it is not deemed a significant detection. GL 688 has only 7 observations, and the periodogram is unable even to identify the correct period.

Seven stars have significant long-term trends, the slopes of which are given in Table 7. These are HR 219a, HR 753, HR 1325, HR 5544a, HR 5544b, HR 7672, and HR 8086. Long-term trends such as these most likely indicate a companion with $`a10\mathrm{AU}`$ and $`M_p\mathrm{sin}i>15M_J`$. Of course many of the other slopes listed in Table 3 may also be due to companions of lower mass (but still with periods $`P11\mathrm{yr}`$). The slopes of HR 5544a and HR 5544b are of almost the same magnitude but opposite sign, as expected for orbiting companions. HR 8086 is an interesting case, because much of its slope is due to secular acceleration, a geometrical effect that stems from the exchange of proper motion for radial velocity. This effect is negligible for the other stars in our sample, but we see it for HR 8086 because it is close and has a large velocity relative to us. Secular acceleration also explains why we don’t see a corresponding negative slope in HR 8085, the companion to HR 8086.

The 20 remaining stars which show significant variability or periodicities are listed in Table 7. These stars are candidates for having planetary-mass companions. We include those stars with a significant periodicity, and those without a significant periodicity but with variability in excess of our prediction of §2.1. We have divided these stars into two groups, chromospherically active ($`P_{\mathrm{rot}}14`$ days) and chromospherically quiet ($`P_{\mathrm{rot}}>14`$ days). The best fit velocity amplitudes for the chromospherically quiet stars are of order the predicted scatter in the velocities due to Doppler errors and intrinsic variability. This makes it difficult, even for very significant periodogram peaks, to confidently identify the observed velocity variations with the Keplerian orbit of a companion. For the chromospherically active stars ($`P_{\mathrm{rot}}14`$ days), there is an additional complication. Even though we account for the excess variability expected in chromospherically active stars, intrinsic velocity variations are likely to be periodic on many different timescales. This renders periodicities seen in these stars suspect, as they may be related to the rotation period, convective motions or magnetic activity. Thus, despite having extremely low false alarm probabilities in the periodogram analysis, none of the periodicities listed in Table 7 are convincing as companions. We are currently making more observations of these stars, to attempt to confirm or rule out these periodicities.

## 4 Upper Limits on Companion Mass

For those stars without a confirmed companion, we would like to know the upper limit on the signal amplitude $`K`$. In this section, we use the periodogram to place upper limits on the velocity amplitude as a function of period. Strictly, our upper limits are on the amplitude of circular orbits because we assume a sinusoidal periodicity. However, we expect our upper limits would not be significantly different for eccentric orbits. This is because for a given period, eccentric orbits have a larger velocity amplitude $`K`$, but an extended duration of roughly constant velocity near apastron. These two effects will cancel each other to some extent for a companion of given mass and orbital radius. Of course, this is not true for orbital periods longer than the duration of the observations, for which our upper limits are strictly for circular orbits. We first describe our method in detail and then present our upper limits on companion mass as a function of orbital radius for each star in the survey.

### 4.1 Method

Our method for placing upper limits uses the fact that, because of measurement errors, different observations of the same signal give different periodogram powers. Most of the time, a large signal amplitude will give a large power. Sometimes, because of noise fluctuations, a large signal will give a small power, but less and less often as the signal amplitude increases. This means that, given a particular observation, we can rule out very large signal amplitudes because they have a very small chance of giving a power as small as the observed value<sup>6</sup><sup>6</sup>6This approach to placing upper limits is a “frequentist” one, as discussed by de Jager (1994) (see also Caso et al. 1998), and applied to the Rayleigh test by Protheroe (1987) and Brazier (1994). It was used in searches for pulsations in the quiescent emission from low mass X-ray binaries by Leahy et al. (1983) and Vaughan et al. (1994). In their work, the noise is dominated by photon counting (Poisson) noise and the time series is evenly sampled. This allows the noise level to be “read off” the power spectrum, giving a natural normalization (Leahy et al. 1983). This is not true in the case of uneven sampling, hence the different approaches in the literature to normalizing the periodogram (see Appendix B).. We define the 99% upper limit to be the signal amplitude such that the periodogram power would be less than or equal to the observed value only 1% of the time.

We proceed as follows. For each star, we find the maximum periodogram power $`z_{\mathrm{max}}`$ from the data (as given in Table 5) for periods between 2 days and 30 years. Then for different trial frequencies, we make simulated data sets of a sinusoid (with frequency $`\omega `$, amplitude $`K`$ and randomly-selected phase $`\varphi `$) and noise. We find the periodogram power $`z(\omega )`$ for each simulated data set and ask how often is $`z`$ larger than the observed value $`z_{\mathrm{max}}`$? The 99% upper limit to the velocity amplitude at frequency $`\omega `$ is that $`K`$ which gives $`z>z_{\mathrm{max}}`$ in 99% of trials. In other words, the observed $`z_{\mathrm{max}}`$ lies at the lowest one percentile of the distribution of $`z`$ that stems from the upper limit to the velocity.

For each trial, we evaluate the periodogram at only one frequency, the trial frequency $`\omega `$. This assumes that the maximum periodogram power will occur at the trial frequency. We have tested this assumption by evaluating the upper limits using 100 frequencies centered on the trial frequency, with the same spacing as used in the periodogram. The upper limits from both methods agree well, and hence we adopt the former as it is computationally faster.

To add noise to the simulated data sets, we utilize the residuals to the sinusoid which best fits the data (i.e. the residuals after subtracting from the data the sinusoid corresponding to the maximum periodogram power $`z_{\mathrm{max}}`$). We motivate this in Appendix B, where we show that, if the data consist of a sinusoid plus noise, the best estimate of the noise variance is the variance of the residuals to the best fit sinusoid. This allows us to obtain the noise variance directly from the data, without having to rely on the estimated errors of §2.1. We find that the variance of the residuals is typically less than the estimated error by a factor of 1.5, consistent with our finding that we overestimated the internal errors (see §2.1). We are thus confident that the residuals give a good estimate of the noise variance for each star.

For each observation, we add noise by selecting at random from the $`v/\sigma `$ distribution of the residuals, scaling by the expected variability $`\sigma _j`$ for each data point. When choosing at random, we select from the pre-fix and post-fix data separately, as appropriate. We have also tried adding Gaussian noise (with appropriate relative weighting between points), scaling by the variance of the residuals to the best fit sinusoid. The upper limits from the two methods agree to 10% or better.

### 4.2 Results

The 99% upper limits are plotted in the $`M_p\mathrm{sin}i`$$`a`$ plane as a solid line in Figure 5 for each star in the sample. For clarity, we state again the meaning of our upper limit. For each star, we imagine repeating the observations many times, each with identical sampling, errors and duration as the real observations. The 99% upper limit is the mass ($`M_p\mathrm{sin}i`$) of a companion which, if present at orbital radius $`a`$, would give a periodogram power larger than that observed in the real data ($`z_{\mathrm{max}}`$) in 99% of our repeated observations. Thus a companion more massive than the upper limit is excluded by the data at better than 99% confidence.

For each star, we calculate the upper limit to the velocity amplitude at 500 logarithmically-spaced periods between 2 days and 30 years. To convert velocity and period $`(K,P)`$ to mass and orbital radius $`(M_p\mathrm{sin}i,a)`$, we use equations (1) and (3). We estimate the mass of each star from its spectral type (Table 2) using a simple empirical formula,

$$M/M_{}=1.3\left(\frac{s}{38}\right)+\frac{3}{10}\left(1\frac{s}{30}\right)^5,$$
(8)   
where $`s`$ parameterizes the spectral type, ranging from $`s=0`$ for F0 to $`s=30`$ for M0 (for example, $`s=22`$ refers to spectral type K2). This formula reproduces Table 9.6 of Lang (1991) to within 10%. Metallicity renders our masses uncertain by a further 10% (see Carney et al. 1994). This is adequate for our purposes, however.

Inspection of the upper limits shows that, for many periods, the upper limit to the velocity amplitude $`K`$ is roughly independent of period. For each star, we list in Table 8 the average velocity upper limit $`\overline{K}`$ calculated for periods less than half the duration of the observations ($`P<T/2`$). We show the corresponding line of constant velocity in Figure 5 as a dotted line. The average velocity upper limit is a good estimate of the upper limit for most periods less than the duration of the observations ($``$ 11 years for most stars, see Table 2).

The upper limit deviates from this “constant velocity” behavior for two reasons. First, the upper limit is larger at periods where the sampling of the observations gives poor phase coverage. For example, many stars show an increase in the upper limit at $`a1\mathrm{AU}`$ because of the tendency for observations to take place at the same time each year. As expected, these aliasing effects are more important in stars with fewer observations. Second, the upper limit to the velocity amplitude increases for periods greater than the duration of the observations. In Figure 5, we show the orbital radius for which $`P=T`$ for each star by a vertical dashed line. The upper limit at long periods ($`PT`$) is sensitive to how we choose the phase in the simulated data sets. The data contain some information about the best fit phase at each period, but we (conservatively) ignore that and choose the phase at random for each trial.

In Figure 5 we do not show results for the 8 stars with confirmed planetary-mass companions (Table 1), or for the 5 stars with companions with $`M_p\mathrm{sin}i>15M_J`$ and $`P20\mathrm{yr}`$ (HR2047, HR5273, HR5553, HR6623 and GL 688). An interesting question is whether these stars have second companions. We are currently generalizing our approach to the case of two companions. For now, we note that the residuals to Keplerian fits to the velocities for the five stars with $`M_p\mathrm{sin}i>15`$ are less than $`30\mathrm{m}\mathrm{s}^1`$ in each case. Searches for second companions to the stars with planetary-mass have been made (see Table 1 for references), without success except HR 3522 ($`\rho ^1`$ 55 Cnc) and HR 458 (Ups And) show evidence for long period second companions ($`P3`$ years).

Walker et al. (1995) (hereafter W95) also used a floating-mean periodogram to place upper limits on companion mass, but took a different approach based on subtracting a sinusoid directly from the data<sup>7</sup><sup>7</sup>7This kind of approach to placing confidence limits is discussed by Lampton, Margon & Bowyer (1976), including the cases where the noise level must be estimated from the data and some of the parameters are “uninteresting” (see also Cline & Lesser 1970, Avni 1976, Cash 1976 and the discussion in Press et al. 1992).. We have applied our method to the data of W95 and find good agreement with their published upper limits. Comparison with their Figure 5 shows that our upper limits on companion mass are about 10–20% lower for $`P<T`$ and 20–30% higher for $`P>T`$. In addition, our upper limits, which show less variability from one period to the next, seem less affected by the sampling of the data. The reason for these differences is not clear. Our comparison shows that there is a “theoretical uncertainty” in the upper limits of about 20%. The good agreement of the two techniques, which are quite different, is encouraging.

## 5 Discussion

So far, we have searched for companions and determined upper limits on the mass of companions for each individual star. We now investigate the implications of our results for the population of planetary-mass companions as a whole.

The observed distribution of the mass and orbital radius of all known planetary mass companions is shown in Figure 9, in which we plot the 17 confirmed companions listed in Table 1 in the $`M_p\mathrm{sin}i`$$`a`$ plane. The dashed lines show constant velocity contours of $`10`$ and $`40\mathrm{m}\mathrm{s}^1`$ for a 1 $`M_{}`$ star. Four features of the observed distribution are interesting. First, all confirmed companions have $`K40\mathrm{m}\mathrm{s}^1`$. Second, no companions have been detected with orbital radius $`a2.5\mathrm{AU}`$. Third, there is a “piling-up” of companions at small orbital radii; for example, of the 17 companions within $`2.5\mathrm{AU}`$, 13 have semi-major axis $`a<0.5\mathrm{AU}`$ and 5 have $`a<0.1\mathrm{AU}`$. Fourth, there is a paucity of companions with orbital radii between $`0.3`$ and $`1\mathrm{AU}`$. In this section, we ask whether our results help to explain these features.

### 5.1 Confirmed and Candidate Companions

First, we summarize the results of our search for companions. We began with a sample of 76 stars from the original Lick Survey. Two stars, HR 5968 ($`\rho `$ CrB) and HR 8729 (51 Peg) were included in the sample because of the discovery of the companions to these stars by other groups (Table 1), and for this reason cannot be considered part of a statistically unbiased sample. Of the remaining 74 stars, six have confirmed planetary-mass companions, or about 8% of our sample. Our periodogram analysis of §3.4 reveals several candidate periodicities (marked “$``$” in Table 5 and listed in Table 7), which are yet to be confirmed or ruled out as being due to companions.

We plot the confirmed companions (circles) and candidate periodicities (squares) in the $`M_p\mathrm{sin}i`$$`a`$ plane in Figure 10. We do this for purely illustrative purposes: we stress that it may be that none of the candidate periodicities are actually due to companions. Open and filled squares indicate candidate periodicities for chromospherically active and quiet stars, respectively. The error bars show the 99% upper limit on $`M_p\mathrm{sin}i`$ obtained as described in §4.1. Broken circles show those confirmed companions not included in our sample of 74 stars. The dashed lines show constant velocity contours of 5, 10, 20 and 40$`\mathrm{m}\mathrm{s}^1`$ for a 1$`M_{}`$ star.

Inspection of Table 7 shows that there are several candidates for companions with velocity amplitude between the Doppler errors, $`K5`$$`10\mathrm{m}\mathrm{s}^1`$, and the lowest velocity amplitude detected, $`K40\mathrm{m}\mathrm{s}^1`$. Thus, one interpretation of our results is that there is an effective detection threshold of $`K40\mathrm{m}\mathrm{s}^1`$. This could be caused by the fact that confirmation of these low amplitude orbits is difficult, as we discussed in §3.5. However, there is an interesting difference in the velocity amplitudes of the candidate periodicities from chromospherically quiet versus chromospherically active stars. Figure 10 shows that there are no candidate periodicities with $`K15\mathrm{m}\mathrm{s}^1`$ in the quiet stars, whereas all of the candidate periodicities in the active stars have $`K15\mathrm{m}\mathrm{s}^1`$. Thus another interpretation is that there is a paucity of companions with velocity amplitudes $`K15`$$`40\mathrm{m}\mathrm{s}^1`$, as none are seen in the subset of quiet stars for which these velocity amplitudes could have been detected. However, we cannot draw conclusions from the present data, as the small number of detections is subject to $`\sqrt{N}`$ fluctuations.

Figure 10 also shows that there are fewer candidates with orbital radii between $`a0.2\mathrm{AU}`$ and $`a1\mathrm{AU}`$ than at other orbital radii. If these periodicities are due to intrinsic stellar variations, a possible explanation is that these variations naturally occur on two timescales, the stellar rotation period ($``$ 1 month) and the timescale of the magnetic cycle ($``$ years). However, the apparent paucity of confirmed companions and candidate periodicities between $`a0.2\mathrm{AU}`$ and $`a1\mathrm{AU}`$ is interesting, and may indicate a real paucity of companions at these orbital radii.

### 5.2 Upper Limits and Detectability

We now turn to the upper limits which we calculated in §4. In Figure 11, we show a histogram (left panel) and cumulative histogram (right panel) of the mean velocity upper limits listed in Table 8. For most stars, we can exclude companions with $`K10`$$`20\mathrm{m}\mathrm{s}^1`$ at the 99% level. Nine stars (15% of the sample) have a mean upper limit $`\overline{K}<10\mathrm{m}\mathrm{s}^1`$ and 38 stars (60% of the sample) have $`\overline{K}<20\mathrm{m}\mathrm{s}^1`$. About 10 stars have $`\overline{K}40\mathrm{m}\mathrm{s}^1`$: inspection of Tables 2 and 8 shows that these stars have either a small number of observations, poor internal errors (for example, because they are faint), short rotation period or a large rms $`>40\mathrm{m}\mathrm{s}^1`$ (for example, because of magnetic activity).

In Figure 12, we plot the cumulative histogram in the $`M_p\mathrm{sin}i`$$`a`$ plane. We take into account the different stellar masses and the effects of sampling (i.e. we use the upper limits as a function of period, Figure 5, rather than just $`\overline{K}`$, Table 8). The solid lines show contours of the number of stars from which a companion of given mass and orbital radius can be excluded at the 99% level (plotted on a 40 x 40 grid). Each contour is labeled by the number of stars from which we can exclude companions above and to the left of the solid line. The dashed lines show constant velocity contours of 5, 10, 20 and 40$`\mathrm{m}\mathrm{s}^1`$ for a 1$`M_{}`$ star. The filled squares show the six confirmed planetary-mass companions in our sample (Table 1). In Figure 13, we show sections of this contour map. We plot the number of stars from which a companion can be excluded (at the 99% level) as a function of orbital radius for different masses, and as a function of mass for different orbital radii.

Figure 13 demonstrates the effect of the finite duration of the observations. Even for massive companions (e.g. $`M_p\mathrm{sin}i>5M_J`$), the number of stars from which a companion can be excluded decreases rapidly for $`a6\mathrm{AU}`$. In Figure 12, we plot with an open circle the point $`M_p\mathrm{sin}i=1M_J`$, $`a=5.2\mathrm{AU}`$, or where our solar system would lie in this diagram if $`\mathrm{sin}i=1`$. This shows that analogs of our solar system are excluded from only a handful of stars once the distribution of $`\mathrm{sin}i`$ is taken into account. This is not true for more massive companions, with $`M_p\mathrm{sin}i>3`$$`4M_J`$. Our results show that these companions are not common at orbital radii $`a1\mathrm{AU}`$, as found by Walker et al. (1995) in their survey of 21 stars. For example, companions with $`M_p\mathrm{sin}i>3M_J`$ ($`>5M_J`$) can be excluded from 80% (95%) of our sample of 63 stars.

Our results give some sense of how the detectability of a companion of a particular mass falls off with radius. For example, for $`M_p\mathrm{sin}i=1M_J`$ ($`0.5M_J`$), the number of stars excluded falls off for $`a0.3`$$`0.5\mathrm{AU}`$ ($`0.1\mathrm{AU}`$). Thus we have detected all companions in the sample with $`a0.3\mathrm{AU}`$ and $`M_p\mathrm{sin}i1M_J`$. For the 40 stars with $`\overline{K}20\mathrm{m}\mathrm{s}^1`$, we can exclude $`1M_J`$ companions out to $`1\mathrm{AU}`$, and $`2M_J`$ companions out to $`4\mathrm{AU}`$. Unfortunately, we cannot assume that the detectability of companions is the same in other surveys. Thus we must restrict our attention to the six companions in our sample.

We have performed a simple test of whether the six companions in our sample could have been drawn from a parent population which is uniformly distributed in orbital radius. First, we assume that the mass distribution of companions is uniform in $`\mathrm{log}M_p\mathrm{sin}i`$ between 0.1 and $`4M_J`$. This seems consistent with the observed distribution of $`M_p\mathrm{sin}i`$ (at least for small orbital radii where detectability is good). We select companions at random from this mass distribution and a uniform distribution in orbital radius from 0 to $`2.5\mathrm{AU}`$. We assign the companion to a star in the sample at random, simulate observations of the star using the observation times and errors from the real data, and ask whether the periodogram power is larger than that observed. The “corrected” orbital radius distribution is then the distribution of orbital radii of those companions which give periodogram powers larger than the observed values. We then use the Kolmogorov-Smirnov (KS) test (Press et al. 1992) to find the probability that the six companions are drawn from this “corrected” distribution.

We find that the probability that the six companions are drawn from a parent population with a uniform distribution in orbital radius is $`10\%`$. If we do not include 70 Vir, which may have a low $`\mathrm{sin}i`$ and thus a mass $`>15M_J`$, the probability is 5%. Thus we cannot rule out that the distribution of companions is uniform in radius. Also, we emphasize again that we are dealing with a small number of detections which are subject to $`\sqrt{N}`$ fluctuations. Future detections will enable us to use the KS test in this way to rule out distributions in mass and orbital radius. For now, without knowledge of the detection thresholds of other surveys, we cannot rule out the possibility that the “pile-up” of companions at low orbital radius or the lack of companions between $`0.2\mathrm{AU}`$ and $`1\mathrm{AU}`$ are due to selection effects or small number statistics.

## 6 Summary and Conclusions

We have presented an analysis of eleven years of precision Doppler velocity measurements of 76 G, K and F type stars from the Lick radial velocity survey. We have performed tests for variability, long term trends, and periodicities. Our sample contains eight confirmed planetary-mass companions, six of which were discovered at Lick, and five companions of stellar or sub-stellar mass (§3.5). Seven stars have significant long term trends, likely indicating a companion with $`a10\mathrm{AU}`$ and $`M_p\mathrm{sin}i>15M_J`$, and 20 stars show variability or periodicities which may indicate a planetary-mass companion (Table 7). We are currently making more observations of these stars, to attempt to confirm or rule out these periodicities. For those stars without a confirmed companion, we have calculated upper limits to the mass of a companion as a function of orbital radius (Figure 5). For most stars, the mean limit on the velocity amplitude is between 10 and $`20\mathrm{m}\mathrm{s}^1`$ or $`0.35`$ and $`0.7M_J(a/\mathrm{AU})^{1/2}`$ (for stellar mass $`MM_{}`$).

We have searched for periodicities and placed upper limits using a “floating-mean” periodogram, in which we fit sinusoids to the data, allowing the zero-point of the sinusoid to “float” during the fit. Allowing the mean to float is crucial to account for statistical fluctuations in the mean of a sampled sinusoid. The traditional Lomb-Scargle periodogram fails when the number of observations is small, the sampling is uneven or the period of the sinusoid is comparable to or greater than the duration of the observations. This may lead to missed detections, or inaccurate upper limits. We have also expanded on the recent discussion by Schwarzenberg-Czerny (1997a, 1997b, 1998) of the correct way to normalize the periodogram and the resulting distribution of periodogram powers. The three different prescriptions in the literature for normalizing the periodogram are statistically equivalent, and all three give a distribution of periodogram powers for Gaussian noise which is significantly different from the usually assumed exponential distribution (see Figure 14). Unfortunately, it is not possible in general to write a simple analytic formula for the false alarm probability, making Monte Carlo methods essential.

Our results help to explain the observed distribution of mass and orbital radius of companions (see Figure 9). The confirmed companions so far have velocity amplitudes $`K>40\mathrm{m}\mathrm{s}^1`$, whereas the Doppler errors lie between $`5`$$`15\mathrm{m}\mathrm{s}^1`$. This is most likely because there is an effective detection threshold which comes about because of the ambiguity introduced by intrinsic velocity flutter (which may be periodic). Confirmation of an orbit is then difficult when the velocity amplitude is similar to the scatter predicted because of intrinsic variability and Doppler errors. We note, however, that in the chromospherically quiet stars ($`P_{\mathrm{rot}}14`$ days), there are no candidate periodicities in the range $`15`$$`40\mathrm{m}\mathrm{s}^1`$. This may reflect a real paucity of companions in this range.

The finite duration of the observations makes it difficult to detect Jupiter-mass companions with orbital radius $`a3\mathrm{AU}`$. Thus the four known companions with $`a>1\mathrm{AU}`$ may be only the first of a population of Jupiter-mass companions at large orbital radii. This is not true for more massive companions, however. It is striking that companions with $`M_p\mathrm{sin}i>3M_J`$ are rare at orbital radii $`4`$$`6\mathrm{AU}`$; we could have detected such objects in $`90\%`$ of stars, yet found none.

A few more years of observations will allow detection of Jupiter-mass companions at $`a5\mathrm{AU}`$, particularly as the poorer quality pre-fix observations become less important. Already, we are able to exclude velocity amplitudes of $`10\mathrm{m}\mathrm{s}^1`$ from 15% of the stars in the sample (for $`a3\mathrm{AU}`$). The velocity amplitudes which can be detected are ultimately limited by intrinsic stellar variability. Even for chromospherically inactive stars ($`P_{\mathrm{rot}}14`$ days), there is intrinsic flutter of a few m/s. Detectability of short period companions should improve in the future as there is more feedback between the observed variability and future observations (ie. stars which show variability are observed more often). Care must be taken to include this effect in assessment of detectabilities.

Our analysis has assumed that the orbits of companions are circular. Yet eccentric orbits appear to be the norm for many of the planetary-mass companions (see Table 1 and Marcy et al. 1999). We find, however, that the periodogram gives a good estimate of the period for eccentric orbits. In addition, we do not expect our upper limits to change substantially for eccentric orbits, as we argued in §4, except for long period orbits for which eccentricities are important. The possibility of a non-zero eccentricity makes identification of an orbital period impossible for periods more than two or three times the duration of the observations. A possible extension of the periodogram to non-circular orbits would be to fit a Kepler orbit at each frequency, and define the periodogram power in terms of the $`\chi ^2`$ of the Kepler fit (see eq. ). It is not clear if the gain in detectability would outweigh the computing power needed for this non-linear least squares fit. A better approach may be to fit higher harmonics as suggested by Schwarzenberg-Czerny (1996).

The observed distribution of companions in mass and orbital radius shows a “piling-up” towards small orbital radii and a paucity of companions between orbital radii $`a0.2\mathrm{AU}`$ and $`a1\mathrm{AU}`$. Because of the small number of companions in our sample, it is not possible for us to say whether these features are due to selection effects. Unfortunately, without knowledge of the detection thresholds of the other radial velocity surveys, we cannot include the other confirmed companions in our analysis. The candidate periodicities we find also show fewer candidates between $`a0.2\mathrm{AU}`$ and $`a1\mathrm{AU}`$, which is intriguing. Future detections will show whether these features reflect the parent population of planetary-mass companions. It may be that companions are to be found at all orbital radii, or it may be that there are two populations of Jupiter-mass companions, one at orbital radii $`a0.2\mathrm{AU}`$ and one at orbital radii $`a1\mathrm{AU}`$.

Either scenario presents an interesting challenge to theorists. Orbital migration models have been proposed to explain the presence of giant planets at small orbital radii (Lin et al. 1999; Murray et al. 1998; Trilling et al. 1998; Ward 1997). These models naturally predict a “piling-up” at small orbital radii (Trilling et al. 1998) because the orbital migration timescale grows progressively shorter as the planet spirals inwards. However, it is not clear exactly how such migration might be halted, in particular at orbital radii as large as $`0.2\mathrm{AU}`$. Indeed, the inevitability of migration may be responsible for the low percentage of solar type stars which have close planetary-mass companions (Ward 1997). If migration depends on gap formation, one would expect migration only to occur for companions above a certain mass. As yet, there is no observed dependence of the mass distribution on orbital radius, except for the lack of companions with $`M_p\mathrm{sin}i3`$$`4M_J`$ at large orbital radius, $`a3`$$`5\mathrm{AU}`$.

There have also been suggestions that gravitational scattering of planets by other planets, companion stars or neighbouring stars in a young star cluster may play a role in determining the final distribution of orbital radii (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Laughlin & Adams 1998; Levison, Lissauer & Duncan 1998). The large range of orbital eccentricities of the observed companions may be evidence for this type of scenario (Marcy et al. 1999). One way to lose enough energy to allow a planet to move from $`5\mathrm{AU}`$ to $`<1\mathrm{AU}`$ may be interaction with the protoplanetary disk during the last stages of dissipation, as suggested by Marcy et al. (1999). It is not known to what extent these different physical mechanisms play a role in determining the distribution of planet masses and orbital radii. Clear theoretical predictions are needed if the discovery of more planetary mass companions is to allow us to distinguish between these different pictures.

We would like to thank Lars Bildsten, Debra Fischer, Evan Scannapieco, Jonathan Tan and Andrew Youdin for useful discussions as this work progressed. Debra Fischer provided rotation periods for several stars, for which a program developed by Phil Shirts was used to measure the Ca IR triplet emission. Jeff Scargle and the referee, Andy Nelson, provided thoughtful comments on the manuscript. Alex Schwarzenberg-Czerny kindly provided a copy of his paper (Schwarzenberg-Czerny 1998). AC was partly supported by the Victor F. Lenzen Memorial Scholarship Fund at UC Berkeley. GWM was partly supported by NSF grant AST 95-20443 and NASA grant NAGW-3182.    

## Appendix A Derivation of the Lomb-Scargle Periodogram

In this section, we sketch the derivation of the traditional Lomb-Scargle formula from a least squares fit of a sinusoid to the data, and show how it relates to the floating mean periodogram used in this paper. Following Lomb (1976), for each trial frequency $`\omega =2\pi /P`$, our model for the velocities is

$$f_j(t_j)=A\mathrm{cos}\omega (t_j\tau )+B\mathrm{sin}\omega (t_j\tau ).$$
(A1)   
The constant $`\tau `$ is introduced to simplify the calculation of $`\chi ^2`$, as the cross terms cancel if we choose $`\tau `$ such that

$$w_j\mathrm{sin}\omega (t_j\tau )\mathrm{cos}\omega (t_j\tau )=0,$$
(A2)   
or equivalently

$$\mathrm{tan}(2\omega \tau )=\frac{w_j\mathrm{sin}2\omega t_j}{w_j\mathrm{cos}2\omega t_j}.$$
(A3)   
The best fit values of the parameters $`A`$ and $`B`$ are those which minimize $`\chi ^2=w_j(v_jf_j)^2`$, where $`w_j1/\sigma _j^2`$ is the weight for data point $`j`$. Setting $`\chi ^2/A=0`$ and $`\chi ^2/B=0`$ and using equation (A2) gives

$`A={\displaystyle \frac{w_jv_j\mathrm{cos}\omega (t_j\tau )}{w_j\mathrm{cos}^2\omega (t_j\tau )}},`$   
$`B={\displaystyle \frac{w_jv_j\mathrm{sin}\omega (t_j\tau )}{w_j\mathrm{sin}^2\omega (t_j\tau )}}.`$ (A4)   
The amplitude $`K`$ and phase $`\varphi `$ of the sinusoid are obtained from $`K^2=A^2+B^2`$ and $`\mathrm{tan}\varphi =A/B`$.

Using equation (A) to substitute the best fit $`A`$ and $`B`$ into $`\chi ^2=w_j(v_jf_j)^2`$, we find the minimum value of $`\chi ^2`$ at each frequency is

$$\chi _{\mathrm{min}}^2(\omega )=w_jv_j^2w_jf_j^2,$$
(A5)   
where $`f_j`$ is now the best fit sinusoid. As in §3.3, we define the unnormalized periodogram power $`\widehat{z}`$ as the reduction in the sum of squares (Lomb 1976),

$$\widehat{z}(\omega )w_jf_j^2=w_jv_j^2\chi _{\mathrm{min}}^2(\omega ),$$
(A6)   
or

$`\widehat{z}(\omega )={\displaystyle \frac{\left[w_jv_j\mathrm{cos}\omega (t_j\tau )\right]^2}{w_j\mathrm{cos}^2\omega (t_j\tau )}}+`$   
$`{\displaystyle \frac{\left[w_jv_j\mathrm{sin}\omega (t_j\tau )\right]^2}{w_j\mathrm{sin}^2\omega (t_j\tau )}}.`$ (A7)   
Equation (A) is the (unnormalized) Lomb-Scargle periodogram (Lomb 1976; Scargle 1982; Horne & Baliunas 1986) modified for unequally-weighted data (Gilliland & Baliunas 1987; Irwin et al. 1989; Scargle 1989). Equation (A6) shows that the floating-mean periodogram we adopt in this paper (§3.3) is a straightforward generalization of the traditional periodogram (see also Walker et al. 1995, eq. \[A2\]).

## Appendix B Normalization of the Periodogram and The Distribution of Noise Powers

There are three different prescriptions in the literature for normalizing the periodogram, namely dividing by (i) the sample variance (Horne & Baliunas 1986; Irwin et al. 1989; Walker et al. 1995), (ii) the variance of the residuals to the best fit sinusoid (Gilliland & Baliunas 1987, this paper), or (iii) the variance of the residuals to the best fit sinusoid at each frequency (Schwarzenberg-Czerny 1996). The motivation for these normalizations is to give the periodogram power $`z`$ a simple statistical distribution when the data are pure noise. The goal is to assess the false alarm probability associated with a given periodogram power. Recent work by Schwarzenberg-Czerny (1997a, 1997b, 1998) has shown that all three normalizations are statistically equivalent, and, at a given frequency, each leads to a simple analytical distribution for Gaussian noise. We now expand on this discussion and describe our choice of normalization.

We point out that because of non-orthogonality between different frequencies, it is difficult to estimate the number of independent frequencies for a given data set. The commonly used empirical formula of Horne & Baliunas (1986) for the number of independent frequencies is based on an inaccurate probability distribution, and is valid only for a particular frequency range. We stress that, unfortunately, it is not possible to write a simple analytic form for the false alarm probability, making Monte Carlo methods essential.

Why normalize the periodogram at all? Scargle (1982) and Horne & Baliunas (1986) showed that if the data points $`X_j`$ are independent Gaussian deviates with variance $`\sigma _0^2`$, the distribution of unnormalized periodogram powers is

$$f(\widehat{z})d\widehat{z}=\frac{1}{\sigma _0^2}\mathrm{exp}(\frac{\widehat{z}}{\sigma _0^2})d\widehat{z}$$
(B1)   
(or simply $`\chi ^2`$ with $`2`$ degrees of freedom). This analysis extends to the weighted periodogram, where $`\sigma _0^2`$ is a measure of the overall normalization of the weights. We would like to stress that if the noise variance were somehow known in advance<sup>8</sup><sup>8</sup>8This is often the situation in X-ray astronomy, for example, where the noise is dominated by Poisson photon statistics, giving a well-defined background power level (Leahy et al. 1983). See also the discussion in Lampton, Margon & Bowyer (1976)., one could simply normalize by the known variance $`\sigma _0^2`$, thus obtaining

$$f(z)dz=\mathrm{exp}(z)dz,$$
(B2)   
an exponential distribution of periodogram powers. However, in many cases $`\sigma _0^2`$ is not known accurately in advance and must be estimated from the data. The idea is to normalize the noise powers to a known level<sup>9</sup><sup>9</sup>9As we noted in §4.1, in principle the background noise level in the power spectrum gives a direct measure of $`\sigma _0^2`$. However, in practice the uneven sampling results in contamination of the noise powers by the signal because of spectral leakage and aliasing effects., aiding identification of localized features in the power spectrum (e.g. periodic signals).

To estimate the noise level from the data set, we use an “analysis of variance” approach (Schwarzenberg-Czerny 1989; Davies 1990). In §A, we showed that

$$w_jv_j^2=\widehat{z}(\omega )+\chi _{\mathrm{min}}^2(\omega ).$$
(B3)   
We now define

$$s^2\frac{1}{Nm}w_jv_j^2,s_f^2\frac{\widehat{z}}{2},s_n^2\frac{1}{Nm2}w_j(v_jf_j)^2,$$
(B4)   
and rewrite equation (B3) in terms of variances, giving

$$(Nm)s^2=2s_f^2+(Nm2)s_n^2.$$
(B5)   
We have partitioned the variance into two pieces, together with their respective degrees of freedom; one piece from the signal $`s_f^2`$ and one piece from the noise $`s_n^2`$. Schwarzenberg-Czerny (1989, 1997a, 1998) showed that under the null hypothesis, $`s_n^2`$ and $`s_f^2`$ are statistically independent by Fisher’s Lemma. This is also true if $`v_j`$ consists of a sinusoidal signal plus Gaussian noise. In this case, $`s_n^2`$ is an unbiased estimate of the noise variance. The different definitions of the normalized periodogram are simply different ratios of the variances in equation (B5). The “traditional” normalization by the sample variance is $`zs_f^2/s^2`$ whereas normalizing by the residuals to the noise gives $`zs_f^2/s_n^2`$.

The distribution of powers when the $`X_j`$ are independent Gaussian deviates can be written down analytically (Schwarzenberg-Czerny 1997a, 1997b). The partitioning of the degrees of freedom means that $`s^2`$, $`s_f^2`$ and $`s_n^2`$ are $`\chi ^2`$ distributed with $`Nm`$, $`2`$ and $`Nm2`$ degrees of freedom respectively. Thus $`z=s_f^2/s_n^2`$ follows an $`F`$ distribution with $`2`$ and $`Nm2`$ degrees of freedom (for example, see Abromowitz & Stegun 1971, §26.6). The distribution of $`z=s_f^2/s^2`$ is more complex as $`s^2`$ is correlated<sup>10</sup><sup>10</sup>10This correlation was neglected by Koen (1990), who incorrectly presumed an $`F`$ distribution in this case. with $`s_f^2`$. The distribution of $`z`$ in this case is an incomplete beta function (see Abromowitz & Stegun 1971, §26.5). The probability that the periodogram power $`z`$ is larger than a given value $`z_0`$ is given in Table 9 for the different normalizations.

Given a probability distribution for the periodogram power, we can write down an expression for the false alarm probability. If the probability that a periodogram power $`z`$ is above some value $`z_0`$ is $`\mathrm{Prob}(z>z_0)`$ (as in Table 9), then the false alarm probability is

$$F=1(1\mathrm{Prob}(z>z_0))^M,$$
(B6)   
where $`M`$ is the number of independent frequencies that were examined. In Figure 14, we show the distribution of maximum periodogram powers for sets of evenly spaced data with $`N=20`$. The crosses are the results of our Monte Carlo simulations. We used three different normalizations for the periodogram. The lines show the theoretical distributions. In the case where the noise variance $`\sigma _0^2`$ is known, we use equation (B6) with $`\mathrm{Prob}(z>z_0)=\mathrm{exp}(z_0)`$ and $`M=N`$ (dotted line). For the normalizations by $`s^2`$ and $`s_n^2`$ (dot-dashed and dashed lines), we use the distributions given in Table 9 and find the best fit value of $`M`$, fitting to the tail ($`\mathrm{Prob}>0.5`$) of the distribution. Both distributions give the same value, $`M=23.5`$. Notice that normalizing by $`s_n^2`$ broadens the distribution of maximum powers (because of the extra uncertainty in the value of $`s_n^2`$), whereas normalizing by $`s^2`$ narrows the distribution (because of the correlation with $`s_f^2`$).

Horne & Baliunas (1986) used Monte Carlo simulations to find the false alarm probability for evenly-spaced data sets<sup>11</sup><sup>11</sup>11Note that due to a typographical error, the values of $`M`$ given in the tables in HB86 are incorrect (Baliunas 1998, private communication).. Normalizing by the sample variance, they assumed $`\mathrm{Prob}(z>z_0)=\mathrm{exp}(z_0)`$ and fit for $`M`$ as a function of $`N`$. However, as Figure 14 shows, the distribution of $`z`$ is different from exponential, especially in the tails of the distribution. Thus the relations of Horne & Baliunas (1986) give inaccurate estimates of false alarm probabilities or detection thresholds. Because the distribution is squashed, the effect is to overestimate both detection thresholds and false alarm probabilities.

In agreement with Press et al. (1992) and Marcy & Benitz (1989), Horne & Baliunas found that $`MN`$ when the period range searched was from the average Nyquist period to the duration of the data set. However, in this paper, we evaluate frequencies several times larger than the average Nyquist frequency. A naive estimate of the number of independent frequencies is $`T\mathrm{\Delta }f`$. Our numerical results show that in general the number of independent frequencies is less than this estimate, hence our Monte Carlo approach for finding the false alarm probabilities in §3.4. It would be useful to have a method, for a given data set, of estimating the number of independent frequencies, allowing one to write down false alarm probabilities analytically. One approach may be to look at the correlations between residuals (Schwarzenberg-Czerny 1991).

Finally, we discuss the difference between normalizing by $`s_n^2`$ at each frequency, and normalizing by $`s_n^2`$ evaluated at the best fit frequency. This choice is really a matter of taste. The distribution of maximum periodogram powers is the same, by definition. However, it seems to us that it is fairer to make comparisons between frequencies using the same normalization for the noise in $`\chi ^2`$. Thus our choice of normalization is that of Gilliland & Baliunas (1987, eq. ): we normalize by the same factor for each frequency, namely $`s_n^2`$ evaluated at the best fit period.