THE FREQUENCY OF LOW-MASS EXOPLANETS. II. THE "PERIOD VALLEY"

The first simulation technique, described fully in O'Toole et al. (2009a), seeks to recover simulated Keplerian signals (with appropriate noise terms arising from measurement uncertainties and stellar jitter) sampled at the epochs of the actual data and retaining the original uncertainties of the actual data. The simulated velocity curves are analyzed using the two-dimensional Lomb–Scargle periodogram (O'Toole et al. 2007, 2009a) which extends the standard Lomb–Scargle periodogram to search eccentricity as well as period space.

The simulated observations for this study were constructed in the following way. We used a grid of periods, eccentricities, and planet masses as input in the same way as O'Toole et al. (2009a), with periods in the range 2–320 days, eccentricities from 0.0 to 0.6 in steps of 0.1, and masses of 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, and 1.0 M_J. Each set of parameters was simulated 100 times, leading to 72,800 simulated velocity curves for each of the 24 stars.

This approach does not use the standard F-test to determine the reliability of a detection. Instead, it uses four detection criteria as follows (O'Toole et al. 2009a). First, the rms of the simulated observations must be greater than or equal to the rms of the residuals of the best-fitting model to those data. That is, the scatter should decrease after the subtraction of a Keplerian model. Second, the measured semi-amplitude must be greater than or equal to twice the measured semi-amplitude uncertainty plus the rms of the residuals of the best-fitting model. Third, the measured period must be greater than or equal to twice the measured period uncertainty. Finally, the value of χ²_ν must be less than or equal to 3. The last criterion is only useful when stellar jitter is not included, as discussed by O'Toole et al. (2009a); in this paper and O'Toole et al. (2009c)—where jitter is included—it is somewhat redundant. The other three criteria are a gauge of the reliability of a fit.

With these criteria, the method does not rely on the assumptions of the F-test. Its chief strengths lie in its ability to recover significantly eccentric orbits (O'Toole et al. 2009a) which would be undetectable using a standard Lomb–Scargle periodogram. In addition, because this method performs a full Keplerian fit to the simulated data, the distribution of orbital parameters derived from the fits provides valuable information on the biases inherent in fitting radial-velocity data (O'Toole et al. 2009a). For the rest of this paper, we will refer to this algorithm as Method 1.

The second technique, described fully in Wittenmyer et al. (2006), adds the Keplerian signal of a simulated planet to the original data and then attempts to recover that signal using a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982). Previous work using similar methodology includes Walker et al. (1995) and Endl et al. (2002). We will refer to this algorithm as Method 2. In this work, we required that the recovered period be within 5% of the injected period, and that its false-alarm probability (FAP) be less than 0.1%. As the null hypothesis of this method is that the data are purely noise and contain no signals due to real planets, the known planets HD 4308b (Udry et al. 2006) and HD 16417b (O'Toole et al. 2009b) must be treated differently. As noted in Wittenmyer et al. (2009), the velocity signature of an unknown planet can be altered by the presence of an additional (known) planet. Hence, we have determined the detection limits for HD 4308 and HD 16417 using the method of Wittenmyer et al. (2009) wherein the simulated planet is added prior to fitting for and removing the known planet's signal. For the Alpha Centauri system (HD 128620/1), we removed a quadratic fit to the data, which provided a satisfactory approximation of the binary orbit over the short time period considered here. We also removed linear trends from the long-period binaries HD 10360 and HD 10361. Each star received seven trials at the 99% recovery level for trial signals with e = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6. In addition, simulations at 90%, 70%, 50%, and 10% recovery were performed at e = 0.0 for each star.

This method is computationally faster than Method 1, since each trial only requires the computation of a simple periodogram, and there are considerably fewer such trials. Method 2 also has the advantage that the noise characteristics of the real data are preserved, whereas for simulated data, additional assumptions must be made about the noise term. However, it breaks down when there are fewer than ∼30 data points, as the reliability of the FAP calculation is strongly dependent on the number of data points. With limited data, even a very strong periodogram peak may fail to reach the 0.1% FAP level. This is less relevant for the current work, as all 24 stars considered here have at least 31 observations (Table 1). We note that Cumming et al. (2008) in their appendix give formulae for the FAP calculation in the regime of small N. For future work, the FAP calculations in Method 1 can be modified in this manner to better handle small N.

For the purpose of this study, we define "detectability" as the fraction of simulated planets which were successfully recovered by the detection criteria employed by each method. To obtain a "big-picture" view of the types of planets detectable with our current data set, we have averaged the results over all stars and all eccentricities. Figure 1 shows the detectabilities from Method 1; the contours indicate the fraction of simulated planets that were recovered. The 99% detection limits obtained in this manner from Method 2 are overplotted as a dashed line. Results for each individual star are shown in Figures 2–13. The detectabilities obtained in this work are highly dependent on the characteristics of the data set; HD 20794 (Figure 5) has a large number of observations with little scatter (N = 78, rms = 3.43 m s⁻¹), whereas HD 28255a (Figure 6) has few data points and considerable scatter (N = 37, rms = 7.5  m s⁻¹). These results highlight the need to perform detection-limit computations on a star-by-star basis (O'Toole et al. 2009a) rather than "whole of sample." Table 3 summarizes the results from both methods for each star. We list the mean velocity amplitude K (averaged over all periods 2–300 days) which was recovered at the 99% and 50% levels for each method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, the results from the two methods are in agreement, though for P < 10 days, Method 1 gives mass limits higher than Method 2 by a factor of 2–4 (cf. Figure 1). The reason for this difference is that Method 1 incorporates the added criterion that the semi-amplitude K of the recovered Keplerian orbit must satisfy K/σ_K > 2. Removing this criterion results in an extremely high proportion of false positives (incorrect detections). Method 2, with a periodogram recovery, has no such constraint on the amplitude, only on the period. Both are more conservative than a human would be in analyzing an individual observed data set by hand. On the other hand, they have a quantitative selection process which means that, given their detection criteria and a set of data, one can use these algorithms to draw robust conclusions about planetary frequency and orbital parameter distribution.

Both Methods 1 and 2 address the question of "What planets could have been detected by these data?" In any such detection-limit determination, the aim is to automate and perform thousands of times what a human investigator would reasonably do to detect a signal (O'Toole et al. 2009a). This process necessarily involves some sacrifices. For example, Method 2 attempts to detect an injected signal using the standard Lomb–Scargle periodogram, an extremely fast computation. However, the Lomb–Scargle periodogram is most efficient at detecting sinusoidal signals; when Keplerian orbits become moderately eccentric, their shapes can become significantly non-sinusoidal (Cumming et al. 2008; O'Toole et al. 2009a; Wittenmyer et al. 2009). To work around this issue, Wittenmyer et al. (2006) restricted the eccentricities of the injected planets to e < 0.6, a range which includes about 90% of known exoplanets. The effects of eccentricity remain evident, however, in the "blind spots" where Method 2 indicates that planets at certain periods are completely undetectable (Figure 14). For each star, there are certain trial periods where the data sampling conspires with an injected eccentric signal, such that an ordinary Lomb–Scargle periodogram does not exhibit a significant peak at that period. Method 2 would then reject the trial signal, even at a large amplitude (K_max = 100 m s⁻¹), and these situations result in allegedly undetectable trial signals. Method 1, however, using the Two-dimensional Keplerian Lomb–Scargle (2DKLS) periodogram (O'Toole et al. 2007) which is much more effective at detecting eccentric signals, would not fall victim to this pathology.

Zoom In Zoom Out Reset image size

The significance of a periodogram peak obtained by Method 2 is computed using the analytic formula given in Horne & Baliunas (1986). However, the FAP obtained by this formula tends to be smaller (i.e., more significant) than that obtained by a bootstrap randomization procedure (Kürster et al. 1997). This is due to the fact that noise in radial-velocity data tends to be non-Gaussian. Thus, a human would visually inspect the periodogram, then perform a bootstrap FAP calculation to obtain an accurate estimate of the believability of any peaks. The bootstrap FAP method, however, requires hours to calculate, rather than the milliseconds it takes for the analytic method, and can therefore not be usefully employed in the analysis of tens of thousands of simulations.

Another necessary concession is the requirement that the recovered period be within 5% of the injected period. The purpose of this constraint is to exclude significant periodogram peaks at aliases of the true period. Spurious periods were also often produced by the lunar cycle (∼30 days) and the duration of the densely sampled "Rocky Planet Search" (48 days). All previously published work using this method (Wittenmyer et al. 2006, 2007, 2009) has used this 5% criterion. This has the effect of preventing incorrect detections (false positives). Method 1 does not have an explicit criterion of this sort, but the four criteria outlined in Section 2.1 serve to keep the rate of false positives below ∼1%–2%; a detailed discussion of false positives is given in O'Toole et al. (2009a).

To summarize, Method 1 seems to be most appropriate for significantly eccentric orbits because the simulated planets are fitted with a full Keplerian model. Since a value of stellar velocity noise (jitter) is added to the simulated data in Method 1 that method is more effective when the radial-velocity noise is well understood. If the stellar jitter is not well understood, it is more appropriate to use Method 2 since the actual velocity data (which include all noise sources) are used.

Jones et al. (2003) and Udry et al. (2003) noted an emerging "period valley" between 10 and 100 days; the ∼250 additional planets discovered since then have continued this trend. Udry et al. (2003) also showed that this feature is mainly attributable to the lack of massive (m sin i > 2 M_Jup) planets in that period range. Our results support this claim, since the 24-star sample contains three low-mass planets with periods between 10 and 100 days (HD 4308b, HD 16417b, and HD 115617c), and the data can exclude planets with m sin i > 0.5 M_Jup at the 99% level for all of this region (Figure 1). The deficit of higher-mass planets with P < 300 days is readily apparent, suggesting that the period valley is not a selection effect; this is consistent with the conclusions reached by Cumming et al. (2008).

At lower masses (10–100 M_⊕), however, the situation is far less clear. This is an extremely interesting mass range, as core-accretion theory, and numerical simulations based on it predict a "planet desert" to arise from rapid gas accretion by cores once they reach about 10 M_⊕ (Ida & Lin 2004, 2008a; Liu et al. 2009). Recent simulations by Liu et al. (2009) showed a pile-up of planets at 0.2–0.3 AU (P∼ 30–50 days). This is a mass range, then, which should show detectable features in the structure of the period distribution. However, this is also a mass range at which serious selection effects impact on the detectability of low-mass planets and most especially in the critical 15–60 day orbital period range. This is why Liu et al. (2009) noted that the predicted pile-up had not yet been observed. Figure 15 shows the period distribution of planets with M_p< 100 M_⊕. A steep falloff in the distribution is evident at P≳ 5 days. This can be understood as arising from two selection effects. First and most obviously, longer-period planets have weaker radial-velocity signals. Second, planet search observing runs are almost exclusively scheduled on large telescopes during bright time, resulting in large peaks in the window function at 28 days. The only means to mitigate this sampling problem is to observe through dark lunations. This is a strategy that only the AAPS is currently employing in its "Rocky Planet Search" observing campaigns and which (as demonstrated by the discovery of 61 Vir c, P = 38 days, K = 3.6 m s⁻¹) has proven successful in finding low-mass planets in this period range.

Zoom In Zoom Out Reset image size

When observing strategies are employed which mitigate selection effects in this period regime, we do indeed find low-mass planets (e.g., Vogt et al. 2010; O'Toole et al. 2009b). These discoveries indicate that the observed falloff in the distribution of such planets (Figure 15) may—contrary to what is seen at higher masses—indeed merely be a result of observational biases. If we look at the results of our simulations, they indicate that 10% of planets with K∼ 2.7 m s⁻¹ could have been detected in our sample of 24 stars. One such planet with an even smaller radial-velocity amplitude has actually been discovered in our sample (61 Vir b, K = 2.2 m s⁻¹), suggesting that at least a further nine more such planets could have been missed and will remain below the detection threshold until additional data are obtained.

In addition, the solar-type dwarf stars typically observed by precision radial-velocity programs have typical rotation periods in the 20–50 day range. Hence, emergent periodicities are not usually believed until sufficient data are in hand to be confident that the periodicity is planetary in origin rather than due to rotational modulation by starspots. This selection bias suggests that low-mass planets with P∼ 30–50 days may be present in extant radial-velocity data sets, but as yet not believed and unpublished.

The importance of planetary migration has been appreciated for some time (Goldreich & Tremaine 1980; Lin & Papaloizou 1986). The modeling work of Trilling et al. (2002) and Armitage et al. (2002) has provided evidence that observational data for exoplanets are consistent with a model where planets are formed at around 5 AU and then undergo migration. This leads to a distribution of planetary semimajor axis a that depends on log(a). Although this scenario provides a good match to observations (e.g., Armitage 2007), migration rates at short orbital periods are expected to be different for different mass planets. The simulations of Rice et al. (2008) indicate three regimes: (1) low-mass planets (<0.5 M_Jup) migrate inward and stall at approximately half the rotation period of the star, (2) high-mass planets (10 M_Jup) should collide with the star, and (3) the intermediate-mass (1 M_Jup) planets population should have an inner edge located at twice the Roche limit (Rasio & Ford 1996). For the purposes of our models, we implement a distribution of planetary semimajor axis that depends on log(a), referred to here as "logarithmic migration."

Although our target star sample and the number of detections are small, we can use the fact that it is well characterized in an attempt to understand the features of the larger inhomogeneous data set of radial-velocity exoplanets. Following the work of O'Toole et al. (2009c) and using our computed detectabilities and a simplistic logarithmic migration (Armitage 2007), Figure 16 shows the expected number of planets for different underlying mass functions as a function of the semimajor axis for the case of equal numbers of exoplanets in each bin. The bold asterisks represent the planets detected from this sample (Table 2). We find that the total expected number of detections for different assumed mass functions observed with our sample varies by a factor of about two for the mass functions shown, with all giving a smooth falloff in the number of detections toward longer periods. Thus, although the peak of the histogram serves to approximately characterize our detection of two exoplanets in the 10–20 day period bin and is consistent with those below the threshold, the clear rise at longer periods in the observed distribution in Figure 16 is not reproduced. However, it must be noted that the sample of all radial-velocity exoplanets shown in Figure 16 is highly inhomogeneous and represents an observed sample of a higher mass range.

Zoom In Zoom Out Reset image size

Although Figure 16 indicates that the distribution in period space is relatively insensitive to different mass functions, we find greater sensitivity when detectability is plotted as a function of mass. Following O'Toole et al. (2009c), we show in Figure 17 the expected detections for our sample as a function of mass for a variety of assumed mass functions. We find that the peak in the expected number of observed exoplanets moves for different mass functions. α = −1.7 leads to a peak in the number of detected planets around 0.02 M_Jup (6.3 M_⊕), the α = −1.2 expected detections peak at around 0.05 M_Jup (16 M_⊕) while the flatter α = −0.7 mass function produces an expected detection peak at around 0.10 M_Jup (32 M_⊕).

Zoom In Zoom Out Reset image size

Given that we have found two exoplanet signals in the 16–25 M_⊕ minimum mass range (corresponding to the 0.05 M_Jup or 16 M_⊕ bin in Figure 17), we may consider these as providing rough limits on the underlying mass function as determined by how changes in the mass function cause the peak of detections to move toward and away from the 0.05 M_Jup (16 M_⊕) bin. We find that the peak moves from the 0.02 to the 0.05 M_Jup bin for α values shallower than α = −1.6 and from 0.05 M_Jup to 0.1 M_Jup at α = −1.1. It is noteworthy that the peak mass for different mass functions is sensitive to migration. If no migration is employed then the mass function peak moves away from the 0.05 M_Jup (16 M_⊕) bin when 0.6 ⩽ α ⩽ 1.2. Based on Poisson-counting statistics of two detections our simulations indicate a normalization of 0.17 ± 0.12 (corresponding to the F value presented in O'Toole et al. 2009c). Our estimate of the underlying mass function shown in Figure 17 is consistent with the first attempt presented in O'Toole et al. (2009c). That is, the additional planet detections in this sample and the comprehensive detection limits presented here serve to strengthen the case for an underlying planet mass function of the form dN/dM ∝ M^−1.0. Given that we robustly do not detect more massive planets in this sample, flatter mass functions are not favored.