TWO PLANETARY COMPANIONS AROUND THE K7 DWARF GJ 221: A HOT SUPER-EARTH AND A CANDIDATE IN THE SUB-SATURN DESERT RANGE^*

After more than a decade of planet search discoveries, the current planetary candidate census has extended to more than 700 planets and counting. These have been discovered using a range of techniques including radial velocity (RV; Mayor & Queloz 1995; Butler et al. 2006; Udry et al. 2007), photometry of transiting planets (Henry et al. 2000; Charbonneau et al. 2000), microlensing (Bennet 2009; Gaudi 2010), direct imaging (Kalas et al. 2008; Maroiset al. 2008; Lafrèniere et al. 2011), and transit timing variations (Nesvorný et al. 2012). To date, no planets have been found using astrometry, although the astrometric signal of already known planets has been detected (Benedict et al. 2002; McArthur et al. 2004; Bean et al. 2007; Reffert & Quirrenbach 2011).

The Kepler mission has recently identified more than 2000 planet candidates (Bathala et al. 2013), including several in the terrestrial mass range. Current and future transiting surveys that have been designed to find habitable planets around nearby stars are the ground-based MEarth project and the Transiting Exoplanet Survey Satellite. Still, ground-based precision RV surveys remain the only technique that has been able to find terrestrial mass planets around nearby stars in their habitable zones (Udry et al. 2007; Vogt et al. 2010; Anglada-Escudé et al. 2012). Recently, much of the RV effort has been directed toward the detection of low-mass planets around the lowest-mass stars: for a given planetary mass and period, a low-mass M-dwarf host will show a larger Doppler amplitude than a G-dwarf host, making lower-mass planets easier to detect.

Currently, the most precise Doppler planet search facilities can achieve an internal long-term precision of 1–3 m s⁻¹ (Mayor et al. 2011; Vogt et al. 2010), which is enough to detect short-period Earth-mass planets around M dwarfs. At this stage, the characterization of the intrinsic-stellar noise sources such as jitter, convective granulation, and asteroseismological p-mode oscillations is becoming more important (Dumusque et al. 2012; Tuomi et al. 2013).

Consequently, the PFS Magellan Planet Search is targeting a sample of 500 of the nearest stars (<100 pc), 200 of which are chromospherically quiet late K and early M dwarfs. The PFS Magellan Planet Search makes use of the Carnegie Planet Finder Spectrograph (PFS). PFS was commissioned at the 6.5 m Magellan Clay telescope at Las Campanas Observatory in Chile beginning in 2009 September and went into full operation in 2010 January. PFS is a temperature controlled high-resolution spectrograph (Crane et al. 2006, 2008, 2010). The spectrograph is maintained at a constant temperature (±0.005°C) near 25°C throughout the year so that the internal optical focus will remain constant without need for adjustment and the refractive index of the air will be stable. An iodine absorption cell (Marcy & Butler 1992) is mounted in front of the instrument's entrance slit, imprinting the reference iodine spectrum directly on the incident starlight, providing both a wavelength scale and a measure of the spectrometer point-spread function (Butler et al. 1996). The iodine cell is a temperature controlled sealed pyrex tube, such that the column density of iodine remains constant indefinitely. PFS was built exclusively to obtain precision Doppler velocity measurements of GKM stars.

With the unexpected diversity of extrasolar planets found to date, every new planet continues to contribute to our understanding of planet formation. In this sense, multiplanet systems are especially important, as they offer the best way to constrain evolutionary history by detailed orbital and dynamical analysis. In this work, we present the detection of two low-mass companions to the nearby low-mass star GJ 221 (BD-061339), by reanalyzing 61 HARPS high-precision RV measurements and including 17 PFS high-precision RV measurements. In Section 2, we describe the observations and data from both spectrometers. Section 3 describes the stellar properties of GJ 221. The orbital analysis is described in Section 4. Finally, a discussion and conclusions are presented in Section 5.

2.1. PFS Observations and Data

Using a 0.5 arcsec slit, PFS obtains spectra with a resolution of R ∼ 80,000 in the iodine region (5000–6300 Å) and covers a wavelength range of 3880–6680 Å. Only the iodine region is used in the Doppler analysis, and the Ca ii H and K lines are used to monitor stellar activity.

Total exposure times range from 300 s on the brightest objects to 720 s on fainter ones. This scheme ensures adequate signal-to-noise ratio (S/N) in the iodine region (S/N >300 per spectral resolution element) and guarantees proper cosmic-ray removal for the faintest targets since longer exposure times yield too many cosmic rays to be handled properly in an automated way. In the case of bright targets, a minimum exposure time of 300 s allows proper time averaging over unresolved low-order stellar p-modes.

Calibrations were acquired at the beginning and at the end of each night. These include 30 flat-field images with the slit illuminated by an incandescent lamp, two incandescent exposures passing through the iodine cell, using both the 03 and 05 slits, and two exposures of a rapidly rotating B star taken through the iodine cell.

Extraction of the spectra from raw CCD images was carried out using an IDL-based reduction pipeline that performs flat fielding, removes cosmic rays, and measures scattered light from inter-order pixels and subtracts it.

The extracted spectra cover the full wavelength range of the spectrograph. No sky subtraction is done as the sky brightness around the iodine region is negligible compared to the brightness of the sources.

Doppler shifts from the spectra are determined with the spectral synthesis technique discussed in greater detail in Butler et al. (1996). In the case of PFS, the iodine region is divided into ∼800 chunks of 2 Å each. Every chunk produces an independent measure of the wavelength, PSF, and Doppler shift. The final measured velocity is the weighted mean of the velocities of the individual chunks. The internal uncertainties presented for all of the PFS RVs are derived as the uncertainty in the mean of the ∼800 velocities from one echelle spectrum. The derived RVs using PFS spectra are presented in Table 1.

Table 1. PFS Velocities for GJ 221

JD	RV	Error
(−2,400,000)	(m s−1)	(m s−1)
55198.67	−1.87	0.93
55581.60	3.77	1.32
55671.48	−6.34	1.61
55786.94	−6.15	1.26
55844.88	4.53	1.27
55845.83	−1.00	1.59
55850.88	0.52	1.45
55851.85	5.68	1.27
55852.86	−0.30	1.88
55853.87	−0.61	1.46
55953.67	−3.59	1.54
55955.63	2.31	1.12
55956.63	0.81	1.24
55957.59	1.57	1.05
55958.64	−4.20	1.58
55959.66	1.13	1.57
55960.66	7.96	1.31
56282.71	−0.19	1.19
56284.69	−5.89	1.13
56291.73	−9.73	1.14
56354.55	7.01	1.34
56356.56	0.00	1.12
56358.55	2.47	1.49

Download table as:  ASCIITypeset image

Since the start of the project, we have monitored a number of stable main-sequence stars with spectral types from late F to early M. Figure 1 shows the RVs of three of these stars during a period spanning more than two years of observations. The Magellan/PFS system achieves measurement precision of 1.5 m s⁻¹, as demonstrated by this figure. The internal median uncertainty of the observations corresponds to the best estimate for the uncertainty due to photon statistics. For the stable stars (rms ∼ 1.5 m s⁻¹ and internal uncertainty of ∼1.1 m s⁻¹), the total uncertainty due to systematic errors, stellar jitter, and unknown planets is therefore ∼1 m s⁻¹ ($\sqrt{1.5^2 - 1.1^2}$). In order to mitigate at least one of these sources of error, stellar jitter, we have selected the most chromospherically quiet stars and have chosen the observing strategy described above.

Zoom In Zoom Out Reset image size

2.2. HARPS Public Spectra

GJ 221 (HIP 27803) has been classified as both a K7 dwarf (Reid et al. 1995) and M0V (Upgren et al. 1972) with an apparent visual magnitude of V = 9.69 and a color B − V = 1.35. The HIPPARCOS parallax (van Leeuwen 2007) gives a distance of 20.31 pc, yielding an absolute magnitude of M_V = 8.15. This makes the star slightly below the zero-age main-sequence relation given in Wright (2004), which would be compatible with a slightly metal-poor main-sequence star (metallicity between −0.5 and the solar one). According to the UVW Galactic velocities (Hawley et al. 1996), GJ 221 should belong to the young disk, giving a lower threshold for its age of 0.1 Gyr. Such a lower limit for the age matches the lower limit obtained from the (V − I_c)–age relations given by Gizis et al. (2002), but it is not a very informative constraint. The chromospheric activity indicator log($R^{\prime }_{\rm HK}$)= −4.8 suggests an age closer to 4 Gyr (Mamajek & Hillenbrand 2008) and the lack of strong flaring events during 70+ observations also suggests an age older than 1 Gyr. As we discuss later, the spectroscopic analysis is in agreement with evolutionary models if an age higher than 1 Gyr is assumed. At least one direct metallicity measurement has been reported before (Casagrande et al. 2008). It reported a metallicity of [Fe/H] = −0.26. The internal scatter on their sample was rather large (∼0.2 dex) so such measurement must be taken with care.

We first estimated the mass of GJ 221 using the mass–luminosity relation by Delfosse et al. (2000) and the V- and K-photometric band values given in Table 3 and obtained M_* = 0.70 M_☉. Assuming this mass, an age of 4 Gyr, and the aforementioned metallicity, we use the evolutionary models in Baraffe et al. (1998) to derive an effective temperature of ∼4300 K. This value is too high compared to previous spectral type determinations (K7–M0 star should correspond to T_eff ∼ 4000 K or lower). In order to double check the stellar parameters with an independent approach, we also derived some of them spectroscopically. The fitting procedure consists of forward modeling synthetic PHOENIX spectra to match the observations using a least-squares solver. The new PHOENIX spectral library has recently been released and is described in full detail in Husser et al. (2013). In addition to flux normalization for the continuum and resolution of the instrument (instrumental line-profile assumed Gaussian), the synthetic spectra also depend on the parameters of interest: T_eff, [Fe/H], log g, and v_rotsin i. The analyzed spectrum was generated by co-adding all the HARPS observations and has an approximate S/N of 400 pixel⁻¹ at 6100 Å. Because atmospheric models of cool stars are still not very accurate in predicting spectra with regions heavy blanketed by molecular features, we focused on small sub-sections of it centered on well-known Fe and Ca i lines (see Figure 2) that are sensitive to log  g and metallicity. The solution was initialized close to our preliminary estimates (discussed before) and the initial value of log g was taken from Takagi et al. (2011, who assumed solar metallicity). The solution converged to log g = 4.5  ±  0.1, T_eff = 4040  ±  50 K, and a metallicity of [Fe/H] =−0.08  ±  0.1. Also, a rotational velocity of v_rotsin i = 1.8 ± 0.1 km s⁻¹ was obtained by fitting the spectra. Using these values, the evolutionary models yield a mass of 0.63 M_☉ if an age of ∼5 Gyr was assumed (any age between 2 and 10 Gyr would fit equally well). This value is ∼11% lower than the value obtained from the mass–luminosity calibration which is known to be uncertain at the 10% level in this mass range (Delfosse et al. 2000). The star is right in the middle of the transition from K to M spectral types and a small change in the parameters leads to a sensible change in derived properties such as luminosities and masses. All things considered, we think that the most consistent picture is recovered when the stellar mass of GJ 221 is obtained from the spectroscopic adjustment + evolutionary models (0.65 M_☉) and that 5% (0.032 M_☉) is a realistic estimate of its uncertainty. The same interpolation of the evolutionary models produces a luminosity of 0.095 L_☉. Updated star parameters and uncertainties are summarized in Table 3.

Zoom In Zoom Out Reset image size

The final orbits were derived with our custom-made software (see description of the methods below) and the solutions were double-checked with the SYSTEMIC package (Meschiari et al. 2009). The detection false alarm probability (FAP) of each candidate was evaluated as follows. First, a least-squares periodogram (Cumming 2004; Anglada-Escudé et al. 2012) was computed on the residuals to the k-planet fit. It is known that when using a periodogram of the residuals, correlations and aliases (e.g., Lovis et al. 2011b) can decrease the apparent significance of additional low-amplitude signals in multi-planetary systems. To account for this, the solution provided by the periodogram of the residuals is refined at each test period by allowing the Keplerian parameters of all of the already detected k planets to adjust. The orbit for the (k + 1) planet candidate is kept circular but its orbital period, phase, and amplitude are also allowed to adjust, ensuring that the optimal multi-planet solution is explored at each test period. The powers (F ratio) of the refined solutions at each test period are shown as red dots in Figure 3. In addition to the planetary signals, the period search model also includes one RV zero point for each data set (γ_HARPS and γ_PFS) and a linear trend (e.g., caused by very long period companions). The F ratio at each proposed period is then computed as

$$\begin{eqnarray} &&F[P] = \frac{\left(\chi ^2_k-\chi ^2_{k+1}[P]\right)/(M_{k+1}-M_{k})}{ \chi ^2_{k+1}/(N_{{\rm obs}}-M_{k+1})}, \end{eqnarray} \tag{ 1 }$$

where M_k is the number of free parameters of the solution including k planets (5 × k from the previous k planets plus one RV offset for each data set, plus one linear trend), and M_{k + 1} is the number of parameters when one additional circular orbit is included (M_{k + 1} = M_k + 2). We call this procedure the recursive periodogram method because of its ability to recursively refine the best Keplerian solution at the period search level. For numerical efficiency purposes, the Keplerian model and the partial derivatives of the observables we use are based on the recipes given in Wright & Howard (2009). Phase folded RV curves of both candidates are shown in Figure 5.

The F ratio of the preferred solution is then used to estimate its FAP analytically using the general methods described in Cumming (2004). In the same study, it was shown that such analytic FAP estimates tend to be overoptimistic because some implicit hypotheses of the method are typically not satisfied by the data (e.g., perfectly known uncertainties, Gaussian posterior distributions, strongly uneven sampling, etc.). This statement was made more precise in Anglada-Escudé & Tuomi (2012), where it was shown that FAP estimates based on the analysis of residual Doppler measurements suffer from strong biases that can work either way (increase false-positive rates, reduce sensitivity to the detection of small signals). Therefore, empirical FAP assessments are always mandatory when periodogram methods are used to assess the significance of a new Doppler candidate. To save unnecessary computations, we will only compute empirical FAPs if the analytic FAP is already found to be lower than 2%. Such empirical FAP computation consists of generating a number of synthetic data sets (or trials) and quantifies how many times a false alarm is induced by an unfortunate arrangement of the noise. Each synthetic data set is obtained by (1) randomly permuting the residuals to the k-planet fit over the same observing epochs (keeping membership to each instrument) and (2) adding back the signal of the k-planet model. We then apply the recursive periodogram method to the synthetic set and save the highest synthetic F ratio in a file. We note that compared to the more classic Lomb–Scarge periodogram of the residuals, each recursive periodogram can be computationally very expensive (for GJ 221, it took 5 minutes to compute a recursive periodogram to search for a second signal). While this is reasonably quick on the real data, it becomes computationally very intensive if 10³–10⁴ trials are necessary to confirm the significance of a detection. Fortunately, each trial is independent of the others, making this task very easy to parallelize. We arbitrarily choose our detection threshold at an FAP of 1%. This threshold keeps the required amount of trials down to a reasonable number while keeping the sample of RV candidates relatively clean of false positives. To save unnecessary computations, we first run 10³ trials. If no FAPs are found, we assume that the FAP is significantly lower than 1% and accept the candidate as significant. If the empirical FAP estimated from the first 10³ trials is higher than 2%, then we also stop the simulations and reject the candidate. Only if the FAP is within 0.1% and 2%, do we then run 10⁴ trials and reevaluate the FAP. If the derived FAP is within 0.5% and 1%, then we extend the simulations to 5 × 10⁴ trials and accept the final result as the empirical FAP. For this system in particular, 10³ trials to search for a second planet running on 40 logical processors were obtained in about 4 hr and 10⁴ were obtained in ∼2 days.

The sampling cadence can also cause confusion between real periodic signals and their corresponding aliases. The modulus of the window function (Dawson & Fabrycky 2010) is shown in Figure 4. The window function is strongly peaked at a few strong sampling frequencies f_s of the time series (related to sidereal day). For a real signal of period P_k, very strong aliases will appear at periods satisfying 1/P_a = |1/P_k + f_s|. The stronger sampling frequencies in this case are f_s =  ±  1.0027, ± 2.0028, ± 3.0055 day⁻¹ (see Figure 4). Note that the first f_s reaches a modulus of 0.81, meaning that a perfectly sinusoidal signal will always be accompanied by two peaks of ∼80% of its power. The exact ratio will depend on interference with additional signals and the spectrum of the noise. The top panel of Figure 3 shows the recursive periodogram search of the combined data sets and the empirical FAP thresholds for the first candidate. Such a periodogram identified a first candidate with an extremely low analytic FAP (<10⁻¹⁰, P ∼ 125 days). As seen in the top panel of Figure 3, it is accompanied by its corresponding aliases around one day corresponding to 1/P_a = 1/125.06  ±  1.0027 day⁻¹. Since no false alarms were found in the first 10³ trials, the candidate was quickly accepted. Using the stellar parameters listed in Table 3, the best Keplerian fit yields a minimum mass of 0.16 M_JUP in a slightly eccentric orbit e = 0.17. The recursive periodogram search for a second planet (Figure 3, bottom panel) revealed a second dominant signal at 3.8741 days. A very prominent alias of this signal is also clearly detected at 1.34285 days (f_s = −1.0027 day⁻¹), 0.79307 days (f_s = +1.0027 day⁻¹), and 0.5731 days (f_s = −2.0028 day⁻¹). The signal is apparent both in the periodogram of the residuals (black line) and in the recursive periodogram refinement (red dots). The corresponding Keplerian solution results in a minimum mass of Msin i of 6.34 M_⊕ in a slightly eccentric orbit (e = 0.17). The empirically estimated FAP required 10⁴ trials and was found to be about ∼0.44%. Therefore, this candidate also satisfied our detection criteria and was added to the solution. The three horizontal lines in both figures represent the 10%, 5%, and 1.0% FAP levels as derived from the empirical FAP calculations. No candidates with an analytic FAP lower than 2% appeared in the recursive periodogram search for a (2+1)-planet solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Confidence intervals for the model parameters were computed using the Bayesian Monte Carlo Markov chain (MCMC) method explained in Ford (2005) but with a slightly different choice of free parameters (see Figure 6). The model parameters describing one Keplerian orbit are always referred to some reference epoch t₀ (which we arbitrarily choose on the first date of the HARPS observations) and they are the period P, semi-amplitude K, orbital eccentricity e, the initial mean longitude λ₀, and the argument of the periastron ω (angle between the periastron of the orbit and the ascending node). The angle λ₀ is defined as λ₀ = M₀ + ω and its choice instead of M₀ (mean anomaly at t₀) is justified as follows. When the orbit is close to circular, the position of the periastron with respect to the plane of the sky is ill-defined and uncertain. As a result, ω and M₀ are strongly degenerate for small eccentricities. Similarly, in the time domain, the precise orbital phase at which the periastron is crossed (instant t for which 2π/P(t − T₀) + M₀ = 2π × N, where N is an integer number) is also very uncertain. As has been discussed in previous works (e.g., Anglada-Escudé et al. 2013), the sum λ₀ contains all the relevant information about the orbital phase when the orbit is close to circular and is much better behaved (non-degenerate) than the two other angles separately. In a technical sense, the convergence of the Markov chains is also greatly improved thanks to the elimination of the degeneracy between M₀ and ω at any e. Let us note that the best-fit values for M₀ can be trivially recovered using ω and λ₀ from Table 4.

We limit our analysis of the activity to searching for signals in the S-index with a period coincident with (or close to) any of the RV candidates. Some other activity indicators are provided by the HARPS-DRS (bisector span, full width at half-maximum of the cross correlation function), but these indices could only be examined on 51 HARPS observations that produced meaningful cross correlation function RV values. In any case, a quick look at their periodograms indicated that none of them showed a periodicity with an analytic FAP lower than 40%.

The HARPS-TERRA software computes the S index in the Mount Wilson system using the prescription given by Lovis et al. (2011a) on the blazed corrected one-dimensional spectrum as provided by the HARPS-Data Reduction Software. In the case of PFS, instrumental S indices were derived from the blaze-corrected spectra following the definition given by Baliunas et al. (1995). We only used the Ca ii H line, because the K line is in a spectral order with very low typical fluxes, which introduces undesirable noise. To convert these values to the Mount Wilson standard system, a simple linear calibration between the PFS and Mt. Wilson systems was made. The main problem of calibrating the S scale is that a number of standard stars observed with both instruments are needed. To address this, we turned to our previously measured values (Arriagada 2011) for some of these stars, performed a linear least-squares fit, and applied the relation to the present measurements.

We computed the least-squares periodogram of the combined series and found a significant periodic signal at 1.84 days (and a corresponding daily alias at 2.16 days) plus strong power at long periods, so we included the adjustment of the trend at the period search level (as we did for the RVs). It is likely that the long-period trend is due to the stellar magnetic cycle of the star (several years?). However, the origin of the 1.84 day signal is less clear. The rotation period of the star could be responsible for such a signal but this would make GJ 221 a very fast rotator and likely to suffer from strong RV jitter (see, e.g., TW Hya in Setiawan et al. 2008). When the 1.84 day signal is removed, a significant periodicity remains at 15 days, again with no RV counterpart. This would be within the range of the expected rotational period of a main-sequence K7 star (Engle & Guinan 2011). In any case, none of the periods detected in the S-index time series appear to have a counterpart in the RVs. Activity induced signals acting at different timescales (1.86 days, 15 days, and several years) are a likely component of the excess noise found in the RVs. In conclusion, although periodic signals are detected in the S index, neither candidate has a period compatible with them, so the Keplerian origin remains the most likely explanation for the signals in the data.

We have analyzed combined precision RV measurements of the K/M dwarf GJ 221 obtained using Magellan/PFS observations and HARPS public spectra. The HARPS-TERRA measurements combined with 17 additional PFS observations reveal the presence of two low-mass companions.

Figure 7 shows our planets in the semimajor axis–mass and semimajor axis–eccentricity parameter spaces of all known extrasolar planets. All of our new detected planets lie well within the parameter space envelope. It is important to note that the longer period planet (GJ 221b) lies in the middle of the so called giant-planet desert predicted by population synthesis models (Ida & Lin 2004, 2010; Mordasini et al. 2012), which is in accordance with other observational evidence (Cumming et al. 2008; Howard et al. 2012). Given the low-metallicity content of GJ 221, the planetary mass of GJ 221c is well in agreement with what observations show in this mass regime (Sousa et al. 2008), as well as what it is predicted by models (Mordasini et al. 2012) in that lower-mass planets do not preferentially form around metal-rich stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Quantifying planet occurrence for M dwarfs provides major constraints for further tunings of these models as well as support for post disk effects such as secular planet migration or planet–planet interactions.

Using the relation provided by Charbonneau et al. (2007), the companion with shorter period (GJ 221b) has a significant probability (8%) of transiting in front of the star and deserves further photometric follow-up. Finally, the star is slightly metal depleted compared to the Sun ([M/H] ∼ −0.26), giving further indication that low metallicity (e.g., Anglada-Escudé et al. 2012) does not inhibit the formation of planets in the super-Earth/Neptune mass regime.

D.M. and P.A. have been supported by the Basal Center for Astrophysics and Associated Technologies (CATA) PFB-06, and by the Ministry for the Economy, Development, and Tourisms Programa Iniciativa Científica Milenio through grant P07-021-F, awarded to The Milky Way Millennium Nucleus. G.A.E. is supported by the German Federal Ministry of Education and Research under 05A11MG3. We are grateful for the advice, support and useful discussions with Ansgar Reiners (IAG), Mathias Zechmeister (AIG), and Hugh Jones (HU). We thank Sandy Keiser for setting up and managing the computing resources available at DTM/CIW. We are grateful for the support that Frank H. Pearl provided at a critical time in the development of the Planet Finding Spectrograph. This work is based on data obtained from the ESO Science Archive Facility under request number GANGLFGGCE175695. This research has made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France, and the NASA Astrophysics Data System.

Note added in proof. After submitting our paper, Curto et al. (2013) submitted an independent detection of these two planetary signals using 102 measurements from the HARPS spectrograph.