© ESO, 2015

1. Introduction

Crossfield et al. (2015) have recently reported from K2 (Howell et al. 2014) observations that the star K2-3 (2MASS 11292037-0127173, EPIC 201367065) hosts a planetary system with three transiting super-Earths. The star is an inactive M0 dwarf (T_eff = 3900 ± 190 K, [Fe/H] =−0.32 ± 0.13 dex), and it is bright enough to be amenable to transit spectroscopy. With a radius of 1.5 R_⊕, Planet d orbits close to the inner edge of the system’s habitable zone. The mass of the planets has not been measured to date.

To constrain those masses, we monitored the radial velocity of K2-3 with the HARPS velocimeter (Mayor et al. 2003). We jointly analysed these new velocities and the K2 photometry through n-body integrations to obtain the mass of planets. Combined with the planet sizes, this constrains their mean densities and bulk compositions and helps explore the rocky-gaseous transition in the super-Earth regime.

2. Observations

2.1. K2 light curve

The K2 mission observed K2-3 during its Campaign 1 (Summer 2014) in long cadence mode¹. We downloaded the pixel data from the Mikulski Archive for Space Telescopes (MAST)² and used a modified version of the CoRoT imagette pipeline (Barros et al. 2014) to extract a light curve. The procedure first determines the circular synthetic aperture that maximises the photometric signal-to-noise ratio on the mean image. For each image, it then computes a modified moment method centroid (Stone 1989) and recentres a heavily oversampled version of the image on the centroid before extracting the flux inside the pre-determined aperture. The degraded pointing stability of the K2 mission couples with pixel sensitivity variations to introduce position-dependent systematics in the raw light curves. To correct for this flux dependence with position, we used a procedure inspired by Vanderburg & Johnson (2014). The light curve was divided into seven equal duration segments, and for each segment we performed a 1D decorrelation as described in Vanderburg (2014). We found that a 21-pixel photometric aperture results in the best photometric precision of the final light curve with a 204 ppm mean rms. The 80-day light curve contains eight transits of Planet b, four of Planet c, and two of Planet d. We only modelled the light curve around those transits, after normalising it with a parabolic fit to its out-of-transit part. To account for the 29.4 min integration time of the long cadence data, we oversampled the model light curve by a factor of 20 and then binned it to the cadence of the data points.

2.2. HARPS radial velocities

We obtained differential radial velocities (RVs) with HARPS, the ESO velocimeter installed at the focus of the 3.6 m telescope at La Silla Observatory (Mayor et al. 2003; Pepe et al. 2003). We chose not to use simultaneous wavelength calibration and to instead rely on the <1  m s^-1 stability of the HARPS spectrograph, since we expected significantly larger photon noise errors. We tried to secure two 1800 s observations per night and collected 66 spectra over a timespan of 103 days.

For an optimal extraction of the velocity signal, we used all recorded spectra (already extracted and wavelength-calibrated by the online pipeline, Lovis & Pepe 2007) to build stellar templates. We shifted each spectrum by its barycentric correction and co-added those shifted spectra into higher S/N templates after carefully masking the telluric absorption lines.

When using such a small set of relatively noisy spectra, we would bias the velocity if we included the spectrum analysed for Doppler shift in the template that it is matched against, since their common noise pattern will contribute to the match. To avoid this bias, we built one template for each spectrum by co-adding all spectra but the spectrum itself.

We then computed the RV shifts that minimise the χ² difference between the individual spectra and their templates (e.g. Zucker & Mazeh 2006; Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015). The 66 resulting velocities, listed in Table A.1, have a 4.9  m s^-1 dispersion for a median uncertainty of 2.9  m s^-1 (compared with 7.1  m s^-1 and 4.3  m s^-1 for the velocities measured by the HARPS pipeline). Although the orbital periods are known from the photometry and do not need to be identified from the RVs alone, Fig. A.1 shows their periodogram for illustration.

2.3. Stellar activity

The HARPS spectral range includes both Hα and the Ca ii H and K lines, which are good tracers of stellar activity. While Hα is a pure absorption line in K2-3 and does not measurably vary, its Ca ii lines, as for all M dwarfs, do have emission cores. We quantify the Ca ii H and K flux through the S index (Vaughan et al. 1978) and find that it varies by 30% on a timescale that is longer than than the ten-day period of Planet b. K2-3 has a lower average S index than Gl 846, an M0.5V star with a ~10.6-day rotation period (Bonfils et al. 2013), consistently with P_rot> 10 days. We used our recent calibration (Astudillo-Defru et al., in prep.) to compute the index from S, and from the vs. P_rot in the same paper, we estimate P_rot ≃ 38 days. The K2 lightcurve varies with a <2 mmag semi-amplitude and a ~40-day characteristic timescale, consistent with a fairly inactive star and with the estimated P_rot. If that photometric variability is entirely due to a dark spot rotating in and out of view, the corresponding RV semi-amplitude is at most ~1.6  m s^-1.

In an effort to quantify activity-induced RV variations, we developed new activity metrics. We built an “active” and a “non-active” template from the third of the observed spectra with the highest and the lowest S-indices, respectively, computed two sets of active and non-active RVs, and used their difference ΔRV = RV_act − RV_noact as an activity tracer. The ΔRV time series (Fig. 1) varies by ±8 m/s for BJD−2 450 000 < 7070, and is consistent with zero after that date. We surmise that the early large variations result from a spot that not only occults part of the star but also imprints the overall spectrum with its spectral signature. That spot would have been present during the first ~days of the observations and then disappeared.

We computed our nominal RVs (Sect. 2.2) with a template that includes all spectra regardless of their S index, and it presumably has a sensitivity to activity that is intermediate between that of the active and non-active template. We account for its stellar activity sensitivity by introducing a term proportional to ΔRV in the RV model, with the proportionality coefficient α adjusted as a free parameter in the fit.

3. Analysis

We analyse this multi-planetary system with an n-body dynamical model that describes the gravitational interactions between all components of the system and not just the pull of the central star. Unlike a Keplerian model, the dynamical model therefore accounts for the photometric and RV changes induced by the mutual attraction of the planets, such as transit timing variations (TTVs), transit duration variations (TDVs), or more generally, transit shape variations. Even though we detect no gravitational interactions, using a dynamical model helps constrain the masses and eccentricities by excluding values for those parameters that would result in detectable interactions or in a highly unstable system.

Our dynamical model uses the well known mercury (Chambers 1999) n-body integrator code to compute the three-dimensional position and velocity of all system components as a function of time. The line-of-sight projection of the stellar velocity is compared to the HARPS RVs. The sky-projection of the planet-star separations are used, together with the planet-to-star-radius ratios and the limb darkening coefficient, to model the K2 light curve (Mandel & Agol 2002). We couple the dynamical model with a Monte Carlo Markov chain (MCMC) code, described in detail in Díaz et al. (2014), to explore the posterior distribution of the parameters. For each MCMC step, we run mercury three times with the same model parameters:

once over the 80 days of theK2 observations, using the Bulirsch-Stoeralgorithm and a 0.02-day timestep to model the light curve with a1 ppm maximum photometric error;

once over the (disjoint) 103 days of the HARPS observations, again with the Bulirsch-Stoer algorithm and a 0.02 days timestep, to model the RVs;

once over a 1000-yr interval, using the hybrid symplectic/Bulirsch-Stoer integrator with a 0.5-day timestep (1/20th of the shortest orbital period).

The last run is not compared to any observation but used to ensure the short-term stability of the system, while rejecting any MCMC step where two orbits intersect or where a planet comes within 0.05 au of the star. (Tidal forces are not included in mercury.) One would ideally like to evaluate stability on longer time scales, with 1000 years chosen as a compromise between ideality and computational expense.

The physical parameters of the model are the stellar mass and radius, the coefficient of a linear limb-darkening law, the systemic velocity, the planetary masses, the planet-to-star radius ratios, and the planetary orbital parameters (a, e, i, ω, n, and M; see Table A.2) at reference time 2 456 812 BJD. To minimise correlations, however, we use a different parametrisation for MCMC jump parameters: orbital period (P), conjunction time of the first transit observed by K2 (T₀), RV semi-amplitude (K), , and . Additionally, we fit a global light curve normalization factor, one multiplicative jitter parameter for each data set and α, the amplitude of the ΔRV term of the RV model (Sect. 2.3).

We use non-informative uniform priors for all MCMC model parameters except the stellar mass and radius, for which we adopt the Gaussian distributions of Crossfield et al. (2015; we also adopt all parameters in common as the starting point of our chains) and the longitude of the ascending node for which we adopt a Gaussian distributions with σ = 2° to enforce the observed physical bias towards coplanarity (Fabrycky et al. 2014). We adopt a non-informative prior for the limb-darkening even though this widens the distributions of several parameters, because limb-darkening models have often been found to be inaccurate when the data quality is high enough to probe them. We ran 48 MCMC chains of 180 000 steps and combined their results as described in Díaz et al. (2014).

4. Results

Table A.2 lists the mode and the 68.3% confidence interval for the system parameters. Figures 1 and A.2 show the RV measurements and the transit light curves, together with their respective models. By fitting a stellar spectral energy distribution to the K2-3 photometry (Crossfield et al. 2015, Table 1) as described in Díaz et al. (2014), we obtain a 47.5 ± 6.0 pc distance and E_{(B − V)}< 0.056 with 68.3% confidence.

Our results agree with Crossfield et al. (2015) for all parameters in common. In addition to adding the RV information, our analysis differs in that we neither impose a circular orbit nor use informative priors on R_⋆/a or the limb darkening coefficient. We also model the transits of all three planets jointly rather than one planet at a time, enforcing consistency in the derived stellar properties. Our approach measures the stellar density more accurately than can be inferred from the average spectrum.

The times of the individual transit are measured with standard errors (which roughly translate to a 1σ detection limits for TTVs) of 1.4, 2.8, and 2.9 min for planets b, c, and d, respectively. Over the short time span of the HARPS observations, the RV of the dynamical model differs from that of a three-Keplerian model by at most 27 cm s^-1, which is well below the ~ 2 m s^-1 measurement noise.

Owing to the low eccentricities and almost edge-on inclinations, all three planets almost certainly undergo secondary eclipses. Table A.2 lists the epochs and durations of these secondary eclipses.

To evaluate the robustness of our mass measurements, we extended the RV by adding a linear drift, a fourth planet, a term proportional to the bisector of the spectra, and any combination of them. We also experimented with restricting the analysis to those RVs that seem unaffected by stellar activity (BJD−2 450 000 > 7070). They always produced similar velocity amplitudes (and masses) for Planet b, but a wide range of values for Planets c and d. We conclude that our data robustly measure the mass of Planet b, 8.4 ± 2.1 M_⊕, but not those of c and d. Many more observations and very careful analysis will be needed to disentangle signatures of Planets c and d from the activity signal. Figure 2 shows the RV signal from Planet b after removing the other terms of the model.

When included in the model, the fourth planet converges to a P_e ~ 100-day period and absorbs part of the residuals seen in Fig. 1 (bottom). Model comparison favours the four-planet model only marginally over the simpler three-planet model, with the Bayesian evidence estimators of Perrakis et al. (2014) and Chib & Jeliazkov (2001) giving odd ratios of 4.6 ± 0.3 and 4.9 ± 2.1, respectively, in favour of the four-planet model. More data will thus be needed to firmly establish whether additional planets orbit K2-3.

Figure 3 adds Planet b to the mass-radius diagram of the known small transiting planets. Planet b must contain rock with at most a 50% water envelope.

Planets around M stars are ideal for characterisation follow-up, thanks to the favourable planet-to-star surface ratio. The other known low-mass planets transiting bright M stars have a gas envelope (GJ436b and GJ3470b, Butler et al. 2004; Bonfils et al. 2012) or may even potentially consist of 100% water (GJ1214b, Charbonneau et al. 2009). Our measurement of the density of Planet b shows that it is either entirely rocky or mostly rocky with a water envelope. It is the first planet in this category, and the brightness of the star makes it a prime target for follow-up.

---

Acknowledgments

We thank the ESO La Silla staff for its continuous support, the HARPS observers (C. Mordasini, J. Martins, and A. Sozzetti, as well as L. Kreidberg for her Mandel & Agol code), C. Damiani for discussions of the dynamics of the system, and T. Fenouillet for his assistance with the LAM computing cluster. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This research made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. N.A.-D. acknowledges support from CONICYT Becas-Chile number 72120460. X.B. and J.M.A. acknowledge funding from the European Research Council under the ERC Grant Agreement No. 337591-ExTrA. NCS acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contracts of reference IF/00169/2012, and POPH/FSE (EC) by FEDER funding through the program “Programa Operacional de Factores de Competitividade – COMPETE”. A.S. is supported by the European Union under a Marie Curie Intra-European Fellowship for Career Development with reference FP7-PEOPLE-2013-IEF, number 627202. S.C.C.B. and O.D. thank the CNES for their grants 98761 and 124378, respectively.

References

Online material

Appendix A: Additional figures and tables

All Tables

All Figures