Radial velocity follow-up of GJ1132 with HARPS. A precise mass for planet 'b' and the discovery of a second planet

GJ1132 is a nearby red dwarf known to host a transiting Earth-size planet. After its initial detection, we pursued an intense follow-up with the HARPS velocimeter. We now confirm the detection of GJ1132b with radial velocities only. We refined its orbital parameters and, in particular, its mass ($m_b = 1.66\pm0.23 M_\oplus$), density ($\rho_b = 6.3\pm1.3$ g.cm$^{-3}$) and eccentricity ($e_b<0.22 $; 95\%). We also detect at least one more planet in the system. GJ1132c is a super-Earth with period $P_c = 8.93\pm0.01$ days and minimum mass $m_c \sin i_c = 2.64\pm0.44~M_\oplus$. Receiving about 1.9 times more flux than Earth in our solar system, its equilibrium temperature is that of a temperate planet ($T_{eq}=230-300$ K for albedos $A=0.75-0.00$) and places GJ1132c near the inner edge of the so-called habitable zone. Despite an a priori favourable orientation for the system, $Spitzer$ observations reject most transit configurations, leaving a posterior probability $<1\%$ that GJ1132c transits. GJ1132(d) is a third signal with period $P_d = 177\pm5$ days attributed to either a planet candidate with minimum mass $m_d \sin i_d = 8.4^{+1.7}_{-2.5}~M_\oplus$ or stellar activity. (abridged)


Introduction
GJ1132 is an M dwarf of our solar neighborhood with a known transiting planet detected by the MEarth survey (Berta-Thompson et al. 2015). Owing to the small size, low mass and low temperature of the parent star (∼ 0.21 R , 0.18 M , 3'300 K), the 1.6-day periodic, 2.6-mmag dips observed in its photometric light curves imply the planet X. Bonfils et al.: Radial velocity follow-up of GJ1132 with HARPS has a size comparable to that of Earth (∼ 1.2 R ⊕ ) and a warm equilibrium temperature (∼ 400-600 K). Being 2-3 mag brighter than most other Earth-size planet hosts seen by, e.g., Kepler, GJ 1132 is an appealing system for followup characterization (Morley et al. 2017;Southworth et al. 2017).
The discovery paper already includes a radial-velocity time-series (RV) collected with the HARPS spectrograph. With an orbital model composed of a single planet and a fixed zero eccentricity we measured an orbital semiamplitude of 2.76±0.92 m s −1 , corresponding to a planetary mass of 1.62 ± 0.55 M ⊕ (Berta- Thompson et al. 2015). Although it favours a rocky composition the constraint is loose given the large mass uncertainty; to the point that even a gaseous composition remains possible in a 3-σ range. In addition to bulk composition, the mass of the planet is also an important parameter to determine the scale height of the atmosphere and constrain transmission spectroscopy observations (Schaefer et al. 2016;Southworth et al. 2017).
If a better mass measurement already makes a strong case for pursuing an intensive radial-velocity follow-up, searching for additional planets is an equally good motivation. With a known transiting planet, we know the system is favourably aligned, and that additional planets have a high chance of being seen in transit as well (Gillon et al. 2011(Gillon et al. , 2017. Actually, GJ1132 was recently used in simulations as an illustrative example of such a strategy (Cloutier et al. 2017b).
This paper reports on our radial-velocity follow-up campaign on GJ1132 with the HARPS spectrograph. We identify a first component with period P d = 177 ± 5 days that we attribute to either an outer planet with mass m d sin i = 8.4 +1.7 −2.5 M ⊕ or stellar activity. After subtracting this first signal, we show GJ1132b is now identified with the sole radial-velocity data. Our time series then reveals another planet with mass m c sin i c = 2.7 ± 0.4 M ⊕ and period P c = 8.92 ± 0.01 days. With an equilibrium temperature of 230-300 K, that super-Earth is located near the habitable zone. Matching our ephemeris with Spitzer observations from Dittmann et al. (2017a), we found that, unfortunately, transits of planet c are largely excluded.

Data
From June 6th, 2015 (BJD=2457180.5) to June 21st, 2017 (BJD=2457925.5), we collected 128 observations with the HARPS spectrograph (Mayor et al. 2003;Pepe et al. 2004), including the 25 measurements already published in Berta- Thompson et al. (2015). We chose the high-resolution mode (R=115'000), used the scientific fiber for the target and the calibration fiber for the sky. In practice, sky brightness is low enough and the second fiber is not used by our pipeline. It only serves as a potential a posteriori diagnostic. Exposure times were fixed to 40 minutes, except for one exposure on June 26th, 2016 (BJD=2457566.5) which was shortened to 1700 sec. Also, note that we discarded a 90th measurement otherwise found in ESO archives. It was just 5 sec and taken on June 11th, 2015 (BJD=2457185.5).
The online pipeline produces extracted spectra calibrated in wavelength (Lovis & Pepe 2007). It also computes radial velocities by doing the cross-correlation with a numerical mask (Baranne et al. 1996;Pepe et al. 2002). We use this initial estimate to shift all spectra to a common reference frame, and we co-add them to build a high signalto-noise reference spectrum. We then refine the radialvelocity determination by finding the best Doppler shift between this reference template and individual spectra (e.g. Howarth et al. 1997;Galland et al. 2006;Anglada-Escudé & Butler 2012;Astudillo-Defru et al. 2015).
The radial-velocity uncertainties are evaluated by measuring the Doppler information content in the χ 2 (RV ) profile, using the formalism of Bouchy et al. (2001) and Boisse et al. (2010). That formalism quantifies the radial-velocity uncertainty by doing a weighted sum over the spectral elements with more weight to the spectral elements with a higher derivative. And since the derivative of a spectrum has a higher variability against noise than the spectrum itself, we instead apply the formula directly to the χ 2 (RV ) profile, which has a few hundred times higher signal-tonoise. For GJ1132, a V=13.5 mag star, we estimate that photon-noise contributes to 2-3 m s −1 to the precision of individual measurements.
In addition to radial velocities, we also measure several activity proxies such as spectroscopic indices (H α , H β , calcium S index and Na; see Astudillo-Defru et al. (2017)). Spectroscopic indices can reveal and trace inhomogeneities at the surface of the star which, animated by the stellar rotation, can contribute to an apparent Doppler shift unrelated to the presence of planets (Bonfils et al. 2007;Robertson et al. 2014).
Our time-series for the RVs, their uncertainties, and for activity proxies is reported in Table 3.

Iterative periodogram analysis
We start with the RV time series (Fig. 1a) and compute its Generalized Lomb-Scargle periodogram (Fig. 1b;Press et al. 1992;Zechmeister & Kürster 2009). We choose a normalization such that a power of 1 at a given period means that a sin-wave fit to the data is a perfect fit (χ 2 = 1) whereas a power of 0 means a sin-wave fit does not improve the χ 2 over that of a fit by a constant. We identify several peaks with a high power: one at the period of the known planet GJ 1132b (∼1.63 day), one near the period 8.9 day, and, the strongest, around 171 days. They have power p=0.25, 0.26 and 0.30, respectively.
To evaluate the significance as a function of the power excess that is measured, we create synthetic data illustrative of time-series with noise only. To preserve the sampling the simulated time-series have the same dates as the original. To preserve the distribution of RVs around their mean, the values are picked by randomly shuffling the original time-series. We compute the periodogram of simulated time-series and measure their power maxima. After many trials we can build a distribution of power maxima that serves as comparison with the power values measured in the original time series. We found that, over 10'000 trials, no simulated time-series has power maxima equal or higher than p=0.30, meaning that the peak seen around P=177 day has a false-alarm probability (FAP) lower than 1/10'000=0.01%; equivalent to > 3.8σ significant.
We note that Berta- Thompson et al. (2015) found a rotation period of 125 days and that a sampling rate of ∼ 1yr could produce an alias at the period 190 day, indistinguishable from the 171-day pic. Nevertheless, the signal appears Left column: RV times-series, before any subtraction (a), after subtracting 1 keplerian fit (d) and after subtraction of a 2-keplerian fit (g). Middle column : Periodograms for each RV time-series on the left column. Right column: RV time-series seen on the left column phase folded to the period of maximum power seen in the periodograms of the middle column. The best sine-fit is superimposed. Dash-otted (resp. dotted) lines in panels b, c, e and h are placed at a power level corresponding to false-alarm probability of 1% (resp. (0.1%).
well sampled when phased with a period of 177 or 190 days. and the 177-day pic is therefore unlikely an alias.
We next reproduce the same periodogram analysis with spectroscopic indices. Figure 2 shows these periodograms for Ca (blue), Na (green), H β (red) and H α (cyan). The period of the stellar rotation is shown with the vertical full line and the periods of the three Doppler signals discussed in this paper are shown with vertical dashed lines. We see broad power excess between 80-300 days with a highest peak at the period of the stellar rotation (∼125 day). Peaks of power are also seen near 175 day, calling for caution on the interpretation of the corresponding Doppler signal. If not due to stellar rotation, magnetic cycle can also induce periodic variations (Gomes da Silva et al. 2011). In Sect. 3.2, we perform a more detailed modelling using Markov Chain Monte Carlo and Gaussian Processes algorithms. By statistically comparing models we show the signal is best describe by a keplerian than correlated noise. Nevertheless, since stellar activity can produce RV variations that match a keplerian (Bonfils et al. 2007), we do not consider this comparison definitely favors the planet interpretation.
We fit and remove a keplerian component to the RVs, and now look at the residuals ( Fig. 1; panel d). The most prominent peak is now that of GJ1132b and RV residuals are well modeled by a keplerian with period P=1.63 day (Fig. 1f). The peak itself has a false-alarm probability of less than 0.01%. The detection is now significant even without prior knowledge on the period from the photometric transits. GJ1132b is therefore confirmed from the sole radial-velocity data.
We continue with a model composed of 2 keplerian orbits. Looking at the residuals we now see strong power excess around P=8.9 days (FAP<0.01%; Fig. 1g), also well modeled with an additional keplerian (Fig. 1i).
The rotation measured for GJ1132 is unambiguously distinct from the last two RV periodicities. This makes stellar activity an unlike culprit of these two short-period signals. We attribute them instead to two orbiting planets, namely GJ1132 b and c.

Joint modeling of planets and correlated 'noise'
Here we apply a second, complementary analysis of the data using a non-parametric Gaussian process (GP) regression model of the correlated RV residuals. GP regression modeling works within a Bayesian framework providing a distribution of functions that model the correlations between adjacent RV measurements following the removal of a mean planetary model containing up to 3 planets in our analysis. This technique has been used recently in the literature to model stellar RV activity thus facilitating the detection and precise characterization of planets around active stars (e.g. Haywood et al. 2014;Rajpaul et al. 2015;Cloutier et al. 2017b;Donati et al. 2016;Bonfils et al. 2018). A complete description of the techniques used to simultaneously model the RV variations with a GP plus keplerian orbital solutions can be found in Cloutier et al. (2017a). Here we briefly summarize the key steps/assumptions used in this work. As mentioned, the broad peak centred around ∼ 175 days in the LS periodogram of the RVs spans the star's photometric rotation period of 125 days (Dittmann et al. 2017a). Therefore the stellar activity might be modulated at approximately the star's rotation timescale. Thus we adopt a quasi-periodic covariance kernel for the GP activity model of the form where t i is the i th BJD in the time-series for i, j = 1, . . . , 128. The GP hyperparameters a, λ, Γ, and P GP describe the amplitude of the correlations, the exponential decay timescale, the coherence scale of the correlations, and the periodic timescale (P GP = P rot in photometry) respectively. Because the GP is intended to model only the stellar activity the posterior probability density functions (PDFs) of the four hyperparameters are trained via Markov Chain Monte-Carlo (MCMC) on a training set which is independent of planetary signals. For this purpose, we used emcee (Foreman-Mackey et al. 2013), a python implementation of the affine-invariant ensemble MCMC sampler (Goodman & Weare 2010). We opt for a training set using the MEarth photometry presented in Berta- Thompson et al. (2015) and used in Cloutier et al. (2017b) to model the stellar RV activity using a GP. Here a well-defined solution is easily found whose periodic term is the photometric rotation period. The resulting marginalized posterior PDFs of the GP hyperparameters λ, Γ, and P GP are then used as priors when next we jointly model the RVs simultaneously with planetary signals and a quasi-periodic GP. The posterior PDFs of these three hyperparameters from training are each approximated by a 1-dimensional kernel density estimation which can then be sampled during the modeling of the RVs. The covariance amplitude a is left as an effectively unconstrained free parameter in the RV modeling. We then model the RVs with one of four potential planetary models. The first model contains the two planets GJ 1132 b and c. The second model contains the same two planets plus the quasi-periodic GP. The third model asummes three planets GJ 1132 b, c, and d. And the fourth model contains the same three planets plus the quasi-periodic GP. In each considered model the orbital period and time of mid-transit of GJ 1132b are assigned Gaussian priors based on the transit results from Berta- Thompson et al. (2015). The orbital period of the putative GJ 1132c (resp. GJ 1132d) is assigned a uniform prior between 8.8 and 9.0 days days (resp. between 120 and 220 days). We adopt uninformative Jeffreys priors between 0-10 m s −1 on each planet's semi-amplitude. This choice of prior was modified and not found to significantly affect the results. The eccentricities e of each planet were sampled indirectly via the jump parameters h = √ e cos ω and k = √ e sin ω where ω is the argument of periastron, a choice that reduces bias toward high eccentricities (Ford 2006).
Next, we conduct a model comparison using timeseries cross-validation (Arlot & Celisse 2010); a computationally less expensive procedure than computing the fully marginalized Bayesian likelihood and one that is independent of the choice of model parameter priors. The unique model parameters for each model-including the GP hyperparameters-are optimized on each of the 107 training sets which contain between 20 and the full dataset size, less 1 (i.e. 127), chronologically-spaced RV measurements. For each split of the data the testing set is the single measurement taken after the final measurement in the training set. The lnlikelihood (ln L) of the testing data given each model optimized on the training set are then computed. Over the 107 splits of the data we compute the median and median absolute deviation per measurement of each model's ln L. Table 1 reports the resulting ln L ratio for various pairs of competing models where each model's ln L is calculated by scaling the median ln L per measurement from cross-validation to the full dataset size of 128 measurements. We find that the 3 planet model is greatly favoured over the 2 planet model with a > 8σ greater ln L. Furthermore the 3 planet + GP model has a marginally better ln L to the 3 planet model (i.e. 0.3σ higher ln L). We therefore conclude that the more simplistic model with 3 keplerian signals and no GP is the model most heavily favoured by the data. Although it favors the planet interpretation for GJ1132(d), we refrain from a definitive con-  Bonfils et al. (2007). The resulting RV model parameters, assuming that the observed RV variations are due to 3 planets plus residual stellar activity which we model with a non-parametric GP, are reported in Table 2. GJ1132 (d) is presented with parenthesis around its d letter to stress it is not accepted as a planet detection.
In the MCMC from which these results were derived we initialized 200 walkers and ran each chain for a duration of approximately 20 autocorrelation times to ensure adequate convergence of the chains. The steps corresponding to the first ∼ 10 autocorrelation times were treated as the burnin phase and discarded. The results are broadly consistent with the results from the iterative periodogram analysis of Sect. 3.1. The resulting marginalized posterior PDFs of the model parameters are shown in Fig. 6.
As a complement, Figures 7 and 8 show phase-folded RVs and residuals with GP regression removed. The uncertainties are not rescaled. After removing the mean GP model, and the best-fit keplerians for GJ1132 b and c, the residual rms is 2.74 m/s.

GJ1132b in context
Compared to the discovery paper, we revise the RV semiamplitude from 2.76 ± 0.92 to 2.85 ± 0.34 m s −1 . This is faster a gain in precision than that expected from the larger number of RVs (0.92/0.48 128/25). We surmise this is because a 3-planet model is a more adequate description of the data. This improves the mass determination by ∼ 45% from 1.62 ± 0.55 M ⊕ to 1.66 ± 0.23 M ⊕ . Together with the radius of GJ1132b measured by Dittmann et al. (2017a), its bulk density then becomes 5.9 ± 1.9 g cm −3 and thus confirms its rocky nature.
The mass-period diagram (Fig. 3) places GJ1132b in context and compares its mass and radius both with other transit detection and with theoretical curves for different bulk compositions (Zeng & Sasselov 2013). The density of GJ1132b appears compatible with a rocky or a denser compositions. With such a diagram, (Rogers 2015) already observed that below ∼1.6 R ⊕ planets are predominantly rocky, i.e. are preferentially found below the mass-radius curves that include significant water or lighter elements in their composition. And indeed, the largest planet that is more than 1-σ away from the rocky curve is Kepler-60b, with a radius of 1.7 R ⊕ (Steffen et al. 2013;Jontof-Hutter et al. 2016).
Conversely, a similar threshold can now be observed with mass: to the exception of 3 planets with very large uncertainties (Kepler -11b, -11f and-177b Lissauer et al. 2011, 2013;Xie 2014), no planet is seen off the rocky curve by more than 1-σ below a mass threshold of ∼3M ⊕ .

GJ1132c, a temperate super-Earth
With a host star's luminosity L = 0.00438 ± 0.00034 L (Dittmann et al. 2017b) and a semi-major a = 0.048 AU, GJ1132c receives about 1.9 times as much flux as Earth from our Sun. Its equilibrium temperature ranges from 232 K for a Bond albedo equals to that of Venus (A=0.75) and up to 328 K for a Bond albedo A=0.
The most recent works that delineate the habitable zone around M dwarfs (e.g. Kopparapu et al. 2016) place the inner edge for GJ1132 around 1.6 times the stellar radiation received by Earth. GJ1132c would thus be considered significantly too irradiated to have liquid water on its surface. Yet the planet remains of considerable interest in the context of habitability. The concept remains poorly understood and will remain so until inhabited worlds are actually found. The habitable-zone's inner edge is thus subject to change with future works. Also, providing future instrumentation were to reach sufficient contrast to resolve the planet from its parent star, probing the existence of an atmosphere would tell us how resilient this atmosphere can be against stellar irradiation, and thus more generally constrain the habitability of M-dwarf planets regardless of the habitability of GJ1132c itself. At a distance of 12 parsecs however (Berta- Thompson et al. 2015), transmission and occultation spectroscopy are probably the sole methods able to resolve such an atmosphere, meaning GJ1132c would be required to transit.

Transit search for planet c
GJ1132c orbits at a distance of about 49±3 R from its host star and, without prior knowledge on the system orientation, the probability to see the planet in transit at inferior conjunction would be ∼1/50. We do nevertheless have prior knowledge on the system orientation since GJ1132 is already known to host a transiting planet with a measured orbital inclination of 88.68 +0.40 −0.33 degrees (Dittmann et al. 2017a). Considering only that nominal value (88.68 degrees), additional planets with strictly co-planar orbits would be seen to transit up to separations of ∼ 43 R , and would be missed beyond. This limit is most probably inside the orbit of GJ1132c and, at first, one may think the prior knowledge we have from GJ1132b may actually nullify the probability of observing any transit for GJ1132c. This is, however, neglecting both the uncertainty on GJ1132c's orbital inclination and possible deviations from perfect coplanarity.
We can instead include both uncertainty and noncoplanarity in our prior. Using the formalism of Beatty & Seager (2010), we make distributions of inclinations centered around 88.68 degrees and with various standard deviations. In Fig. 5, one can see the prior probability that GJ1132c undergoes transit quickly jumps to ∼43% for a distribution of inclinations with only a small standard deviation of 1 degree.
From our analysis in Sect. 3, we derived an ephemeris for the passage of GJ1132c at inferior conjunction (BJD=2457506.02±0.34). As seen in Fig. 4, this falls inside the long 100-h monitoring made by Dittmann et al. (2017a)  with Spitzer, which covers epochs between BJD=2457502.5 and 2457506.8 with almost no interruption and a sensitivity down to planets smaller than Mars (Dittmann et al. 2017a). Possible transits of GJ1132c are largely ruled out with < 1% chance an existing transit was missed.
The prior probability joined to the probability left by the incomplete coverage of the transit time-window gives

Stellar activity or a cold planet
Based on activity diagnostics, we were not able to rule out stellar activity as the main cause of the 170-d radial velocity periodicity. However, if later observations were to confirm the planetary nature of that signal, GJ1132 d would have a minimum mass m d sin i = 8.4 +1.7 −2.5 M ⊕ , in a mass domain between super-Earths and mini-Neptunes. Its semi-major axis a d = 0.35 ± 0.01 AU would place this planet beyond the ice line, with an equilibrium temperature of 86 K (resp 111 K) for a Bond albedo of 0.75 (resp. 0.3). Planet d would be ∼7.3 times further away from the star than planet c. A priori, its transit probability would thus be ∼ 7.3 smaller. The Spitzer photometry rejects only a small fraction (13%) of possible transit configurations. 1 0.43% is the upper limit considering a standard deviation of 1 degree in our above calculations and for all other standard deviations the posterior probability is less than 0.43%.

Transit-timing variations
The orbital periods ratio of planets b and c is close to 11/2. Even if inside a resonance, its low order would imply low transit-timing variations. We use the Rebound code (Rein & Liu 2011) with the WHFast integrator (Rein & Tamayo 2015) to compute TTVs for this system. Although not shown here, TTVs were found to be generally less than 30s, in agreement with the lack of TTVs found in Dittmann et al. (2017a).

Conclusion
To conclude, our HARPS radial-velocity follow-up helps to bring a more complete, more complex description on the GJ1132 system. We confirmed the detection of planet 'b' with the sole radial velocities, we refined its characteristics including its mass, density and eccentricity. We also detect at least one new planet in the system. Assuming coplanarity with planet 'b', GJ1132 c is a super-Earth with mass m c = 2.75 +0.76 −0.61 M ⊕ . Its equilibrium temperature is that of a temperate planet although it is probably too close to the star to allow for liquid water on its surface. Finally, we also detect a third keplerian signal but leave its true nature -planet or stellar activity-undecided. If confirmed as a planet, GJ1132 d would be a super-Earth or mini-Neptune found further away from the star and beyond the ice line. Further observations, and potentially at other wavelengths with infrared spectrographs such as Carmenes (Quirrenbach et al. 2014;Sarkis et al. 2018) or SPIRou (Delfosse et al. 2013), may help further interpret the true nature of that signal.