THREE PLANETS ORBITING WOLF 1061

The NASA Kepler mission has demonstrated that a high fraction of low-mass stars are planet hosts (Dressing & Charbonneau 2013, 2015; Kopparapu 2013), and the evidence is mounting that multi-planet systems are common with discoveries like GL 581 (Udry et al. 2007), GL 667C (Delfosse et al. 2013), and GL 876 (Rivera et al. 2005). Multiple long-term programs have demonstrated the utility of Doppler techniques for exoplanet detection and characterization (e.g., the HARPS, CORALIE, and ELODIE groups; the Anglo-Australian Planet Search group, the Keck HiRes and Lick Observatory groups). Several new instruments are being built to expand the use of Doppler techniques, with a focus on M dwarf planet hosts (e.g., CARMENES, Alonso-Floriano et al. 2015; SUBARU HDS, Snellen et al. 2015; ESPRESSO Pepe et al. 2014; Veloce).¹

M dwarfs are good Doppler targets due to their many sharp molecular absorption features and their typically slow rotation speeds. Low-mass stars are also low-luminosity stars, which contracts their habitable zones to short periods (typically <100 days). These close-in planets have a stronger Doppler effect on their host star and that Doppler amplitude is further increased by the low mass of the host star. All of these characteristics make M dwarfs excellent targets for rocky, habitable-zone exoplanet searches, e.g., Berta-Thompson et al. (2015).

The 2017 launch of the NASA Transiting Exoplanet Survey Satellite mission (Ricker et al. 2015) will deliver a quantum leap in the number of M dwarf planet candidates for Doppler follow-up—many of them in the habitable zones of their hosts.

The HARPS spectrograph was purpose-built for planet hunting (Mayor et al. 2003). It is one of the most precise exoplanet hunting facilities with a demonstrated long-term precision of <1 m s⁻¹ (Pepe et al. 2003; Dumusque et al. 2015). The extensive HARPS database contains long-term monitoring data for more than 100 M dwarfs and is an excellent resource for the refinement of planet search techniques.

As noted in Dumusque et al. (2015), the HARPS DRS velocities suffer from a significant yearly variation due to the barycentric movement of the spectra across pixel discontinuities in the HARPS detector. A quick analysis of the HARPS DRS data shows that this annual signal is clearly the largest variation present for Wolf 1061.

We have developed a process to extract Doppler velocities from the DRS spectra with improved precision over the standard analysis using a cross-correlation template generated from the data itself, rather than from ab initio line lists. The template generated is therefore inherently differential. We obtained a typical velocity precision of <1 m s⁻¹. Most interestingly, analysis using this template delivers velocities that do not show any significant yearly or half yearly periodicities for Wolf 1061.

The observations and methods used to determine Doppler velocities of Wolf 1061 are outlined in Sections 2 and 3, while Section 4 discusses the velocity and stellar activity results and detected planetary signals. Section 5 focuses on details of the star's habitable zone, and Section 6 provides a brief conclusion of our main results.

We have used the 148 publicly available HARPS spectra of Wolf 1061. The spectra span 380–680 nm at resolution ∼110,000. The observations were taken over 10.3 years and typically have 900 s exposure times and a median signal-to-noise ratio (S/N) of 65 pixel⁻¹ at 6000 Å.

In the most recent HARPS M dwarf sample publication (Bonfils et al. 2013), Wolf 1061, though analyzed with a much smaller number of observations than presented here, was not considered exceptional. No significant periodicities were detected, though a test for variability indicated possible variation at 67.3 days, which we identify as a likely planet orbit.

Details on the calibration and precision from the HARPS DRS reduced data are discussed in Pepe et al. (2003). We use the archival spectra of Wolf 1061 in their 72 order, flux-per-camera-pixel form, i.e., not in their merged, flattened, and rebinned form. Our method of extracting velocities is similar to the HARPS DRS except for three changes:

• 1.  
  an additional cleaning of outlying pixels in the spectra,
• 2.  
  a custom template built for the star,
• 3.  
  an iterative solution using the measured velocities to improve the custom template.

As in the standard HARPS analysis (Pepe et al. 2002), we determine a velocity from each echelle-order spectrum in each observation by cross-correlation against a weighted spectral mask. However, rather than using an ab initio line list, we construct a custom mask for each target star using a high S/N spectrum, which we obtain from the data for that star.

First, all spectra are shifted to the barycentric reference frame and then any known Doppler velocity is removed. In the first iteration the Doppler velocity is unknown and so initially the HARPS DRS velocities are used. In later iterations the computed velocity is used instead. Errors in the spectra are estimated as the square-root of the flux at each pixel. All spectra are then flattened by the spectrograph blaze function (defined as part of the DRS reduction) and rebinned onto a standard axis of 0.01 Å bins. The spectra are then inspected to identify and correct outliers at the 5σ level using the ensemble of data at each wavelength.

Next a very high S/N spectrum is obtained from summing all the rebinned spectra. To construct the mask, this summed spectrum is upsampled to 0.001 Å and inverted. The inversion is achieved by assuming a flat pseudo-continuum at the maximum value in the summed spectrum and subtracting the spectrum from the pseudo-continuum value. This allows precise identification of the depth (now seen as line "height" due to the inversion) and position of all of the pseudo-line absorption profiles present in the summed spectrum. A cross-correlation template can now be constructed from the list of positions and depths so identified.

To obtain Doppler velocities the spectra go through a similar process to the above (i.e., shifting and cleaning), though without the subtraction of the Doppler velocity. The spectra are then rebinned to a logarithmic wavelength scale in 300 m s⁻¹ bins. The spectral template for each order is built using delta functions with the positions and depths obtained above and covers the spectrum up to 2 Å from the ends of each order. There are four wavelength regions considered too contaminated by telluric absorption to be useable for this method: 5850–6030 Å; 6270–6340 Å; 6450–6610 Å; 6855–6910 Å. Finally, the template and spectra are cross-correlated and the cross-correlation profiles for each order are added. The S/N in the summed spectrum is insufficient to construct a reliable template for Wolf 1061 at the blue wavelengths, so the bluest 25 orders are not included in the combined cross-correlation profile.

Velocities are computed by fitting a Gaussian to the combined cross-correlation profile. The velocities are then fed back into the building process for the cross-correlation template as mentioned above. Iterations continue until the difference between velocities from iterations is less than 20 cm s⁻¹.

We compute our Doppler velocity errors in a similar manner to the DRS system and so obtain similar error estimates, typically <1 m s⁻¹. This error is based on the slope and flux present in the combined cross-correlation profile and calculated in the same way as outlined by Butler et al. (1996). Figure 1 shows all of the resulting Doppler velocities.

Zoom In Zoom Out Reset image size

The Doppler velocities obtained from the above method do not suffer the annual variation observed in the DRS velocities. By building the template from the shifted and combined spectra the impact of the "seams" in the CCD (Dumusque et al. 2015) are built into the template. The template we construct is an inherently differential measurement and therefore is less sensitive to the causes of the yearly variation in the DRS velocities, which cross-correlates spectral data with an imperfect wavelength solution against an ab initio line list with a "perfect" wavelength solution. A subsequent paper will discuss our method in more detail to investigate this in greater depth.

We begin with a sequential extraction of periodicities using the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) on the full velocity data set. The significant signals found are used as start points for a genetic algorithm analysis of the velocity data. The top panel of Figure 2 shows the periodogram of the full data set and the location of the three significant periods detected. At each stage, we fit a Keplerian to the data with a starting period corresponding to the strongest periodogram peak. After the subtraction of these three signals no significant variation is observed in the residuals, which have an rms of 1.86 m s⁻¹. We determined the false-alarm probability (FAP) of each signal using a bootstrap randomization approach (Kürster et al. 1997) coupled with the error-weighted generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009). We used 10,000 bootstrap realizations to derive the following FAPs for the three signals: 17.87 days – <0.0001, 67.27 days – <0.0001, 4.89 days – 0.0052.

Zoom In Zoom Out Reset image size

A variety of phenomena including stellar activity and stellar pulsation can induce signals in the RVs (e.g., Queloz et al. 2001; Robertson et al. 2014). Many early M dwarfs have long rotation periods (>20 days) and also have activity, e.g., starspots or magnetic effects, that cause variations in the stellar spectrum that are detected as Doppler velocity changes and hence the rotation period or its harmonic may be misinterpreted as a planet (McQuillan 2013). To identify the detected velocity signals in Wolf 1061 as planets we must rule out the possibility that they are caused by intrinsic variability in the host star or that they are associated with the rotation period of the star.

By following Suárez Mascareño et al. (2015) we can determine the log₁₀(R^'_HK) index for Wolf 1061 from the Ca ii H&K lines which then provides an estimate of the rotation period. Suárez Mascareño et al. (2015) define the so-called S-index using a 0.4 Å window around 3933.664 and 3968.470 Å to define the Ca ii H&K emission and a 20 Å window around 3901.070 and 4001.070 Å for the continuum definition. The computation of log₁₀(R^'_HK) from the S-index was then straightforward using B − V = 1.566 for Wolf 1061 from Landolt (2009). We obtain log₁₀(R^'_HK) = −5.45  ±  0.13 indicating Wolf 1061 is very inactive. Suárez Mascareño et al. (2015) provide an equation for the relationship between the rotation period and the log₁₀(R^'_HK) value but it is poorly constrained for M dwarf stars. Although the relationship provides a value for the rotation period of 76 days, due to its large relatively uncertainty we take from this calculation only that the measured log₁₀(R^'_HK) is indicative of a large relatively long rotation period, e.g., >60 days.

Additionally there were 537 epochs (8.6 years) of precision (error < 4 mmag) V band photometry available from the All-Sky Automated Survey² (Pojmanski & Maciejewski 2004) photometry. Examination of the periodogram of the photometry revealed no significant periodicities (FAP < 0.1) at any period. This is further evidence in support of a rotational period longer than our detected RV periods for Wolf 1061.

Recent work by Rajpaul et al. (2016) has also shown that small peaks in the window function due to the data sampling can masquerade as significant signals after the subtraction of any detected activity signals. Essentially, the subtraction of signals from a data set can suppress other peaks in a power spectrum and leave behind spurious signals. Therefore we examine the following activity periodograms in Figure 2 directly without subtraction.

In addition to the S-index outlined above, there are several other indices known to vary with stellar activity or pulsation: the bisector span, the cross-correlation function FWHM, Hα index, and Na D line index (Vaughan et al. 1978; Santos et al. 2002; Hatzes et al. 2010; Suárez Mascareño et al. 2015). Any significant variation observed in these indices at periods coinciding with detected signals in our velocities would suggest that the signal is likely due to stellar activity, pulsation, or rotation, and not an exoplanet.

We have used a typical bisector span measurement i.e., the difference between the average of the 10%–40% and 60%–90% depths of the bisectors (Pepe et al. 2003). The Hα index emission region was difficult to define because the Hα line is weak and there was very little variation observed in the spectra. We settled on using a 2 Å window centered at 6562.810 Å for the Hα line and a 4 Å window centered on 6556.000 Å for the continuum. We note that the so-called "continuum" as used here is not pure continuum emission (pure continuum is not easily accessible in M Dwarf stars) but rather just a region of stable spectrum that should not vary under the same circumstances that the Hα line varies, e.g., due to starspots or magnetic effects. In this way we measure the strength of the Hα emission and its changes relative to a stable region of spectrum. For the NaD index we used a 0.4 Å window around the D lines at 5895.924 and 5889.954 Å, and a 14 Å window centered on 5875.000 Å for the continuum region. The FWHM of the cross-correlation profile and its error are a direct result of the Doppler velocity measurement process.

In Figure 2 the window function for the periodograms shows no peaks at our velocity signal periods but does have two small peaks at 75 days and 92 days and shows that our long period sampling >300 days is very poor. In Figure 2 panels 3–7 we show periodograms of the activity indices with the positions of the three signals found in our velocity data overplotted.

There are no indications of significant signals in any of the indices at the detected velocity periods. In fact only the S-index shows any significant signals at periods <300 days. The S-index shows several significant peaks at periods >75 days. We interpret the long period peaks (>300 days) in all of the activity indices as likely to be due to the data sampling, since the window function clearly indicates that we are not sensitive to such periods. The several significant S-index peaks at >75 days may be associated with the rotation period of Wolf 1061 but it is not possible to identify what that period might be. Although there is no direct link between this S-index variability and our RV signal periods, it is impossible to completely rule out a connection between our longest RV period (67.2 days) and the S-index signals since rotation related velocity signals could show up as harmonics of the rotation period, i.e., P_rot/2 or ${P}_{{\rm{rot}}}$/3.

We also check for correlations between the velocities and activity indices and obtain correlation coefficients ranging from −0.11 to +0.12. We use a simple two-tailed t-test and find no significant correlations between the velocities and any of the indices examined here (significant p-values are p < 0.05).

Based on these tests Wolf 1061 is a very stable star and hence the signals detected in the RVs are best explained as planets. Taking the RV data set periodicities as starting points, we explored a broad parameter space using a genetic algorithm to fit three Keplerian orbits to the data (e.g., Horner et al. 2012; Wittenmyer et al. 2013). Approximately 10⁷ three-planet models were explored, and convergence occurred rapidly, favoring a system with three low-mass, low-eccentricity planets. We derived final system parameters using Systemic 2 (Meschiari et al. 2009) and estimated parameter uncertainties using the bootstrap routine therein on 10,000 synthetic data set realizations. The parameters given in Table 1 represent the median of the posterior distribution and 68.7% confidence intervals. Using the host star mass of 0.25 M_⊙ (Maldonado et al. 2015), we derive planetary minimum masses of 4.25 ± 0.37 M_⊕  (planet c), 5.21 ± 0.68 M_⊕  (planet d), and 1.36 ± 0.23 M_⊕  (planet b). Figure 3 shows the velocity time-series phased to the period of each individual signal with the other planets removed.

Zoom In Zoom Out Reset image size

Table 1.  Wolf 1061 Planetary System Parameters

Parameter	Eccentricity Free			Circular Orbits
	Wolf 1061b	Wolf 1061c	Wolf 1061d	Wolf 1061b	Wolf 1061c	Wolf 1061d
Period (days)	4.8876 ± 0.0014	17.867 ± 0.011	67.27 ± 0.12	4.8871 ± 0.0011	17.870 ± 0.005	67.28 ± 0.08
m sin i (M⊕)	1.36 ± 0.23	4.25 ± 0.37	5.21 ± 0.68	1.33 ± 0.21	4.10 ± 0.31	4.97 ± 0.58
K (m s−1)	1.29 ± 0.22	2.70 ± 0.25	2.24 ± 0.30	1.27 ± 0.20	2.53 ± 0.19	1.97 ± 0.23
Eccentricity	0.0 (fixed)	0.19 ± 0.13	0.32 ± 0.16	0.0 (fixed)	0.0 (fixed)	0.0 (fixed)
ω (degrees)	⋯	37 ± 47	102 ± 70	⋯	⋯	⋯
Mean anomaly (degrees)	61 ± 41	313 ± 59	170 ± 41	279 ± 30	191 ± 20	255 ± 21
a (AU)	0.035509 ± 0.000007	0.08427 ± 0.00004	0.2039 ± 0.0002	0.035507 ± 0.000006	0.08428 ± 0.00002	0.2040 ± 0.0002
${\chi }_{\nu }^{2}$	3.72	⋯	⋯	3.77	⋯	⋯
rms of fit	1.86 m s−1	⋯	⋯	1.87 m s−1	⋯	⋯
FAP of signal	0.0052	<0.0001	<0.0001	⋯	⋯	⋯
Transit probability	14.0%	5.9%	2.6%	⋯	⋯	⋯
Transit depth (mmag)	2.6	3.3	5.2	⋯	⋯	⋯

Download table as:  ASCIITypeset image

We have detected three signals in the velocity time series at P = 17.867 days, P = 67.274 days, and P = 4.888 days. The variability observed in the activity S-index is at longer periods than our detected velocity signals and does not show any clear association with them, though we acknowledge a possible association with our longest period velocity signal of P = 67.274 days cannot be completely ruled out—for the present we cautiously interpret this signal as an exoplanet.

Our velocities and activity indices and the HARPS DRS velocities are available in Table 2.

Wolf 1061 is a bright and close M dwarf with V = 10.1, $d=4.29\;\mathrm{pc}$, and spectral type M3/3.5 V (Bonfils et al. 2013; Maldonado et al. 2015). Bonfils et al. (2013) provided several fundamental parameters for Wolf 1061, though a recent study by Maldonado et al. (2015) has reported it as being slightly cooler than Bonfils et al. (2013). The stellar parameters taken for the computation of the mass of the planets and the habitable zone of the star are those of Maldonado et al. (2015): effective temperature ${T}_{{\rm{eff}}}$ = 3393 K, mass M = 0.25 M_⊙, and luminosity L = 0.007870 L_⊙.

We use the habitable zone calculations of Kopparapu et al. (2013a, 2014)³ , which give both conservative and optimistic habitable zones (depending on the approach accepted). We obtain conservative habitable zone limits of 0.092–0.18 AU and optimistic limits of 0.073–0.19 AU. These boundaries place the middle planet well inside the optimistic habitable zone, and just outside the inner boundary of the conservative habitable zone. The outer planet in the system lies just outside the outer boundary of the optimistic habitable zone for the host star.

We have found strong Doppler signals in data for Wolf 1061 that indicate the presence of three potentially rocky planets: a 1.36 M_⊕ minimum-mass planet with an orbital period of 4.888 days (Wolf 1061b), a 4.25 M_⊕ minimum-mass planet with orbital period P = 17.867 days (Wolf 1061c), and a probable 5.21 M_⊕ minimum-mass planet with orbital period P = 67.274 days (Wolf 1061d). With such short-period planets, we can consider the possibility that one or more may transit. The a priori transit probabilities are as follows: planet b—14.0%, planet c—5.9%, planet d—2.6%. If the planets are transiting, the planet masses are equal to the minimum masses above. We can then estimate planetary radii from the mass–radius relation of Weiss & Marcy (2014). This yields the following rough radius estimates: 1.44 R_⊕ (b), 1.64 R_⊕ (c), and 2.04 R_⊕ (d). Using the stellar radius given in Maldonado et al. (2015) (R = 0.26 R_⊙) we estimate transit depths of 2.6, 3.3, and 5.2 mmag for planets b, c, and d, respectively. Such signals can be detected from ground-based telescopes in excellent conditions, and we encourage purpose-built facilities like MEarth and MINERVA (Nutzman & Charbonneau 2008; Swift et al. 2015) to pursue the Wolf 1061 transit windows when this target becomes observable again in early 2016.

The 17.867 day planet is of particular interest because it is of sufficiently small mass to be rocky and is in the habitable zone of the host M dwarf star. The probable outer planet at 67.274 days resides just on the outer boundary of the habitable zone and is also likely rocky. These planets join the small but growing ranks of potentially habitable rocky worlds orbiting nearby M dwarf stars.

Based on data obtained from the ESO Science Archive Facility under request: Duncan Wright 189972. This research has been supported by ARC Super Science Fellowships FS100100046 and ARC Discovery grant DP130102695. We also wish to acknowledge the many helpful comments from the referee in producing this manuscript.