EDP Sciences
Free Access
Issue
A&A
Volume 608, December 2017
Article Number A63
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201731307
Published online 07 December 2017

© ESO, 2017

1. Introduction

In spectroscopic time series data, we can measure the radial velocity (RV) variations of a star as a result of the Doppler shift. Orbiting planets induce periodic signals originating from the wobble of the star around their common barycenter. Other sources of RV variations are related to different kinds of surface phenomena induced by short-term effects like stellar oscillations and granulation (Meunier et al. 2015; Dumusque et al. 2011) and the magnetic activity of the host star. This includes mid-term effects like dark spots and bright faculae moving with the stellar rotation (Haywood et al. 2016; Herrero et al. 2016; Suárez Mascareño et al. 2017c) and includes long-term effects connected with the location or number of surface active regions and the magnetic cycle of the star (Lanza et al. 2016; Santos et al. 2010). The changes in RV result most importantly from the different amount of photons coming from the blueshifted and from the redshifted limbs of the rotating star, and also from the effects of convective blueshift and its suppression in magnetically active regions (Haywood et al. 2016; Livingston 1982).

In recent years exoplanet surveys have paid considerable attention to M dwarfs since they provide the possibility to detect, characterize, and understand Earth-like planets. For Sun-like stars, the current technical RV limit of ~1 m s-1 is still an order of magnitude too large to accomplish this goal. In addition to the advantage of their larger RV signals, M dwarfs have closer habitable zones with shorter orbital periods and higher transit probabilities. However, different surface phenomena rotating approximately at the stellar rotation period induce signals into the RV data that can mimic planets that are seen in stars like GJ 581 (Joiner et al. 2014; Robertson et al. 2014), GJ 667C (Anglada-Escudé et al. 2013), or GJ 674 (Bonfils et al. 2007), but can also hide planets of Earth-mass (Howard et al. 2014).

The exoplanetary community is investing significant effort into disentangling magnetic activity and planetary signals and into understanding and correcting for the influences of magnetic phenomena on the RVs. These include computational efforts to improve the search for periodicities in RV data (see, e.g., Dumusque 2016; Dumusque et al. 2017, for an overview), the construction of instruments observing at near-infrared wavelengths like CARMENES (Quirrenbach et al. 2014), GIARPS (Claudi et al. 2016), SPIRou (Artigau et al. 2014), or HPF (Mahadevan et al. 2012), and performing spectroscopic surveys tailored to M dwarfs (Bonfils et al. 2013; Tuomi et al. 2013; Irwin et al. 2015; Alonso-Floriano et al. 2015; Anglada-Escudé et al. 2016).

With our HArps-n red Dwarf Exoplanet Survey (HADES) our goal is to detect rocky planets by monitoring 78 bright M0 to M3 stars. This is a coordinated effort of Spanish (IEEC-CSIC, IAC) and Italian (INAF) institutions. The observations began in August 2012 and, as of March 2017, we have obtained approximately 3700 spectra. We have already detected two low-mass planets around GJ 3998 (Affer et al. 2016) and a potentially habitable one around GJ 625 (Suárez Mascareño et al. 2017b). We have also described in detail the stellar rotation characteristics of our sample by analyzing their activity indices (Maldonado et al. 2017; Scandariato et al. 2017). Additionally, we performed simulations to predict the outcome of our survey and found an average underlying noncorrelated activity jitter of 2.3 m s-1 (Perger et al. 2017).

In this work we search for planetary companions around GJ 3942 (HIP 79126) by analyzing its RV time series measurements obtained with the HARPS-N spectrograph. GJ 3942 is a low-mass star with spectral type M0, located in the constellation of Draco and with kinematics suggesting membership to the young disk population (Maldonado et al. 2017). The analysis of its spectroscopic properties (Maldonado et al. 2017; Reiners et al. 2012; Morales et al. 2008) reveals solar metallicity, intermediate rotation, and average to weak magnetic activity. Just recently the star was labeled as a new candidate RV-variable star based on seven RV data points in the near-infrared obtained by Gagné et al. (2016). We provide an overview of the basic properties of GJ 3942 in Table 1.

Table 1

Stellar parameters of GJ 3942.

In Sect. 2 we present the spectroscopic and photometric observations obtained and give a first overview of the periodic signals present in the data. In Sect. 3 we explain the search for planets in detail and use various state-of-the-art approaches to correct for the effects of magnetic activity. This leads to the detection of a super-Earth sized planet GJ 3942 b and a second planetary candidate that we discuss in Sect. 4. The conclusions of our work are given in Sect. 5.

2. Data reduction and preliminary analysis

We obtained optical spectra with the Northern High Accuracy Radial velocity Planet Searcher (HARPS-N, Cosentino et al. 2012), connected by fibers to the Nasmyth B focus through a Front End Unit of the 3.58 m Telescopio Nazionale Galileo (TNG) in La Palma, Spain. It is a fiber-fed, cross-dispersed echelle spectrograph with a spectral resolution of 115 000, covering a wavelength range from 3830 to 6900 Å. We observed with fixed integration times of 900 s to obtain data of sufficient signal-to-noise ratio (S/N> 20) and to average out potential short-term periodic oscillations of the star (Dumusque et al. 2011), although they do not seem to be present in M dwarfs (Berdiñas et al. 2016).

We reduced the raw data with the YABI tool (Hunter et al. 2012; Borsa et al. 2015), which implements the DRS data reduction pipeline. It uses the classical optimal extraction method by Horne (1986) and includes bias and background subtraction and flat fielding to deliver cosmic ray-corrected, wavelength-calibrated spectra. Furthermore, it calculates the cross-correlation function (CCF), which is the correlation of a spectrum with an M2-type template mask in velocity space. Various parameters, including the RV, are calculated from this CCF. For more details see Cosentino et al. (2012) or the DRS manual1. In our case, we decided not to use RVs determined by the DRS, but instead use the Java-based Template-Enhanced Radial velocity Re-analysis Application (TERRA, Anglada-Escudé & Butler 2012). It handles the full process of unpacking the HARPS-N archive files in the DRS output. The template-matching algorithm has been shown to deliver more accurate RVs in the case of M-type stars (Perger et al. 2017).

GJ 3942 was observed on 145 occasions resulting in spectra with an average S/N of 56.2. The full time series dataset spanning 1203 nights or 3.3 yr is provided in Table 9. As a quality check, we searched for correlations between the RVs and the S/N (which would be indicative of charge-transfer inefficiency effects), for obvious flare events (seen as emission in lines that have a contribution from the stellar chromosphere), or other irregularities and did a 5σ clipping on the RV measurements. As a result, we rejected three data points with S/N< 20 from further analysis. With the 142 remaining spectra, we determine an RV scatter and mean error of 6.01 and 1.13 m s-1, respectively. For comparison, these values are 6.56 and 1.92 m s-1 in the case of the DRS. We also obtain an absolute average RV value of 18.7047 ± 0.0066 km s-1.

2.1. Instrumental radial velocity drift

thumbnail Fig. 1

All available 8199 instrumental drifts distributed over 356 nights. The solid line indicates the mean value and the dashed lines the 1σ interval.

Open with DEXTER

The HARPS-N instrument implements the possibility of using calibration light simultaneously with the target to measure the instrumental drift of the RV during science observations. A second fiber can collect light from a stable hollow-cathode lamp of ThAr, resulting in a simultaneous wavelength solution. Until August 2015, we did not use this calibration strategy because for texp> 200 s it was supposed to contaminate the RV measurements from the science fiber. As we learned during our program, this does not seem to be the case, at least for the brighter objects (V< 10.5 mag) like GJ 3942. Therefore, of the 142 measurements, only 41 have simultaneous calibration from which instrumental drift values can be determined. However, measurements taken by the Italian Global Architecture of Planetary Systems (GAPS, Covino et al. 2013; Desidera et al. 2013) project included simultaneous calibration, and this provided the opportunity of estimating instrumental drift corrections for a total of 51 additional nights on which GJ 3942 was observed.

To estimate drift corrections we made use of 8 199 measurements available from Spanish and Italian observations. They are distributed over 356 nights with an average of 23.03 measurements per night and are shown in Fig. 1. In some nights, the RV drift varied by as much as 12 m s-1. We wanted to use the data of a specific night to estimate the instrumental drift at the time of the observation of GJ 3942. To evaluate the best method, we randomly took out one measured value and calculated the predicted value from the remaining measurements from the given night using several approaches. We show the mean values and rms of the differences in Table 2 using 100 trials. For the 51 nights covered by instrumental drift values, we find that the best results are obtained by averaging the available measurements in a time interval of ±1 h. Thus, we applied this methodology (method 4 in Table 2) and added quadratically an uncertainty of 0.466 m s-1 to each measurement. For those 50 measurements where no RV drift calculation was possible we considered an additional uncertainty of 1.025 m s-1 added in quadrature, which corresponds to the total rms of all drift values. The instrumental drifts for each observation are provided in Table 9.

Table 2

Different methods used to estimate instrumental drifts and the resulting rms of the data.

2.2. Activity indices

In order to distinguish signals induced by activity from genuine planetary Doppler signals, we calculated activity indices of magnetically sensitive features, which are primarily associated with emission from the stellar chromosphere. These are the central cores of the Ca ii H & K lines and the Hα line (Gomes da Silva et al. 2011).

The S index measures the relative flux of the Ca ii H & K emission lines compared to a local continuum. The lines are formed in the hot plasma of the chromospheres of stars, hence vary with the strength of the stellar magnetic field (Anglada-Escudé & Butler 2012), and measure the lower chromosphere (Gomes da Silva et al. 2011). We calculated the S index with the definitions by Duncan et al. (1991). This consists of adding the total flux in the central cores of the two lines and calculating the ratio to the sum of fluxes corresponding to two triangular windows on either side of the lines: S = α·(FH + FK) / (FR + FV). We further modified the wavelength windows following Gomes da Silva et al. (2011) to avoid adding information from the undesired photospheric contribution in the wings of the line, slightly shifted the continuum windows, and used a value α of 18.4 following Lovis et al. (2011). This method produces values that are not in the same scale as the Mount Wilson index (Duncan et al. 1991; Suárez Mascareño et al. 2015), but we are more interested in defining a precise indicator that is useful to search for small periodic variations. The Hα index is calculated as Hα = FHα/ (FR + FV) (Boisse et al. 2009), modified by using a broader central window to include a greater contribution from the chromosphere following Gomes da Silva et al. (2011), and slightly different continuum windows.

Uncertainties for the indices are calculated as Gaussian errors using Poisson noise on the spectral data. The S and Hα indices show relative errors of 1.8% and 1.6%, respectively. The index calculated from the Na i doublet (Díaz et al. 2007) was not used in the present study because our preliminary analysis did not show it to be a good activity proxy for this target.

The DRS pipeline and the CCF technique also provided us with additional quantities. Magnetic phenomena on the stellar surface affect the shape of the lines and this is imprinted on the CCF. Variations of the parameters that describe the CCF may therefore be correlated with activity. Three of these parameters are calculated automatically by the pipeline: the bisector inverse slope (BIS, Gray 1989), the full width at half maximum (FWHM), and the peak value (contrast) of the Gaussian fit to the CCF. The shape of the CCF of our target shows two prominent side lobes caused by the high density of spectroscopic features and is not well fitted by a Gaussian (e.g., Affer et al. 2016; Suárez Mascareño et al. 2017a). This explains the higher dispersion of the DRS RV values and we therefore prefer not to use the CCF indices as activity proxies in our analysis.

2.3. Signal identification

The classical approach to finding periodic signals in unevenly sampled time series data is the Lomb-Scargle (Lomb 1976; Scargle 1982) or the Generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009), which is a sinusoidal fit on error-weighted, shifted values using the function (1)The parameters c0, c1, and c2 are fitted to the data for each period (P) step, and a time shift t0. In the case of a Keplerian orbit, this fit can be specified by (2)where γ is the RV offset; K the semi-amplitude of the signal; e and ω the eccentricity and argument of periastron of the planetary orbit, respectively; and TA the true anomaly of the planet (Wright & Howard 2009).

As a measure of the significance of a possible signal, the false alarm probability (FAP) is calculated by bootstrapping data with 10 000 permutations. Common practice indicates that a signal at the 1% FAP level is defined as suggestive, while a 0.1% FAP level signal is considered to be statistically significant. An analytical formula for the FAP is given by Horne & Baliunas (1986). We analyzed the RVs, the S index, the Hα index, and their respective residuals using the GLS.

thumbnail Fig. 2

Time series data of (from top to bottom) radial velocities (in m s-1), S index, and Hα index (in 10-2) are shown in panels a, c, and e, respectively, as observed (black dots) and their linear trend is indicated by the blue lines. The mean errors of each dataset are illustrated by the red line on the lower left side of each panel. Panels b, d, and f show the de-trended data (i.e., linear fit subtracted) folded to the most prominent periodicity at 16.3 days and with zero phase corresponding to the first observation (black dots) and a simple sinusoidal fit (blue curve).

Open with DEXTER

The left panels of Fig. 2 show the different measurements corresponding to the RVs, S index, and Hα index, where the four observing seasons are evident from the clustering. The global trends of the data yield negative slopes with modest significance. For the S and Hα indices, slopes of 0.039 and 0.040 per year are found, respectively, while the overall slope of the RVs is insignificant at 0.12 m s-1 per year. We subtract these trends from the data for subsequent analysis to correct for any long-term effects resulting from a possible magnetic cycle (Dumusque et al. 2017). We note that a more detailed treatment of the trends leads to similar results. In the right panels of Fig. 2 we show the different de-trended datasets phase folded to the most prominent periodic signal at 16.27 days (the mean of the best fits for the three values) and the best sinusoidal fits. The origin of phases is set to day one of our observations in all three cases. We observe a phase shift of approximately 120 deg between the activity indices and the RVs. This is not unusual and, e.g., Santos et al. (2014) found a shift of 15 deg for their Ca ii index possibly connected to the location of activity phenomena on the stellar surface.

thumbnail Fig. 3

Generalized Lomb-Scargle periodograms of the de-trended time series data and of the residuals after correcting for the 16.3-day signal with a simple sinusoidal fit. In panels a and b this is shown for the radial velocities, in panels c and d for the S index, and in panels e and f for the Hα index. Additionally, we show in panel g the window function of the time series. Yellow lines indicate the important periodic signals for GJ 3942, namely 16.3, 10.4, and 6.9 days, and we also show additional periodicities mentioned in the text. The dashed blue horizontal line shows the 1% FAP level of the respective dataset.

Open with DEXTER

In Fig. 3 we show the GLS periodograms of the three de-trended time series and their residuals from frequencies of 0.0 to 0.3 day-1 including all peaks clearly above noise level. In panel (a) we show the periodogram of the RVs with significant signals at 16.3 and 6.9 days. Together with the signals at around 10 days, these are the main periodicities in our data and they are marked by the yellow bands in the different panels. We correct for the most prominent signal in the RV data at 16.29 ± 0.22 days by a simple sinusoid (see Fig. 2), analyze the residuals with the GLS, and show the periodogram in panel (b). The subtraction cleans up the periodogram significantly, but preserves the 6.9 day peak, with a FAP below 0.1%, as well as the two different peaks at 10.1 and 10.4 days, which increase in significance to reach 1% FAP. Additionally, at the same FAP level, a signal at 5.4 days emerges. The S index periodogram in panel (c) does not show any significant signal besides that at 16.3 days. Only some excess power is observed at approximately 9.9 and 42 days, which disappear after removing the period of 16.26 ± 0.22 days, as seen in the periodogram in panel (d). There, signals at 100, 7.5, and 9 days increase in power, with the latter two likely related to the first harmonic of 16.3 days. The strong signal at 16.3 days is also found in the Hα index dataset in panel (e). Here, a suggestive periodicity of approximately 200 days is present as well. Removing the signal at 16.25 ± 0.22 days leaves only a peak connected to the first harmonic of the 16.3-day periodicity in the periodogram in panel (f). The window function (WF) shown in panel (g) does not contribute to any of our important signals.

This most prominent periodicity in our RV data at 16.3 days is clearly related to rotational modulation and may have suffered variations over the observational time span (Santos et al. 2014; Dumusque et al. 2017). To investigate this, we performed sinusoidal fits to the RV data, but grouped them into seasons. We considered seasons 1 and 2 (S12) together to increase the number of measurements in such a way that 43, 53, and 46 data points were used for S12 (green dots), S3 (red dots), and S4 (black squares), respectively. Fits with a constant period of 16.29 days are shown in panel (a) of Fig. 4. In the diagram, we can observe a phase shift of S3 of approximately 35 deg, while S4 shows an increase in the semi-amplitude of approximately 70%. Fitting a free period also shows its variability since for S3 alone the main signal is at a lower value of 15.9 days. The same trend is also seen for the activity indices in panels (c) and (d).

thumbnail Fig. 4

Illustration of changing characteristics of signals in RV, RV residual, and activity index data of GJ 3942. It is grouped into seasons 1 & 2 (S12, green dots), season 3 (S3, red dots), and season 4 (S4, black squares). The best purely sinusoidal fit for each dataset is shown. Panel a shows the RV data folded in phase with a 16.29-day period; panel b the residuals after correcting for the best sinusoidal 16.3-day period fit and phase folded with 6.91 days; panel c the S index data folded in phase with a 16.26-day period; and panel d the Hα index data phase folded with 16.25 days.

Open with DEXTER

The variation of the prominent signal at least on a yearly basis over the observed time span is another strong argument for the signal being induced by activity phenomena present on the stellar surface and rotating with the star (Lanza et al. 2016; Díaz et al. 2016; Meunier & Lagrange 2013). The changing characteristics of the fit are likely explained by the changing locations, sizes, and strengths of the magnetically active regions and by differential rotation. The periodogram peak at 6.9 days is not affected by the subtraction of the 16.3-day signal and does not show significant power in any of the activity indices. We therefore conclude that this RV periodic signal must be a result of the Keplerian motion of an orbiting planet, GJ 3942 b. In panel (b) of Fig. 4, where the RV residuals are shown, no phase shift in S3 can be seen, but an increase in the signal’s amplitude is visible, whereas the parameters in general seem to be more stable than for the 16.3-day period. While this planetary signal is quite robust, the interpretation of the signals between 7.5 and 10.5 days is more complicated. In this interval, the periodicity of approximately 8.2 days corresponds to the first harmonic of the activity signal (Boisse et al. 2011). The peaks around 10.4 days, which are not highly significant, do not have a straightforward interpretation and the different possibilities are discussed in Sect. 3.

thumbnail Fig. 5

FAP development for the 16.3-day (black) and the 6.9-day (red) signals of the RV data using an increasing number of data points. The dashed horizontal line shows the analytical 0.1% FAP level, and the two vertical dash-dotted lines divide the data into the different seasons.

Open with DEXTER

In Fig. 5 we show the time evolution of the FAPs of the 16.3- and 6.9-day signals in the RV data. We calculate the analytical FAP level of the most prominent peaks inside an interval of 16.3 ± 1.5 and 6.9 ± 1.0 days, respectively, using an increasing number of data points. The 16.3-day signal shows a steady increase in significance over the entire observation time. The planetary signal reached a 1% FAP level with approximately 60 observations, was hidden behind the noise until approximately 90 observations, and finally fell rapidly below a 0.1% FAP level with 100 observations. The diagram shows changes in the significance of the signals between seasons 3 and 4, in agreement with the change in the activity pattern discussed above.

thumbnail Fig. 6

Photometry of GJ 3942 including normalized data from APACHE (panel a), magnitude differences (GJ 3942 and calibration stars) of EXORAP filters B (panel b), V (panel c), R (panel d), and I (panel e), and visual magnitudes from the WASP (panel f) and Hipparcos catalogs (panel g). The red lines on the lower left show the rms of the nightly averages (panels a, f, g) and the RV rms of the calibration stars (b to e) as proxies for the uncertainties of the measurements. We note that APACHE was used during seasons 1 and 2 of the HARPS-N observations, EXORAP from season 2 to 4 and beyond, and WASP and Hipparcos data are from around 2007 and 1991, respectively.

Open with DEXTER

2.4. Photometry

As support to the HARPS-N spectroscopy, GJ 3942 was monitored photometrically by the programs A PAthway toward the Characterization of Habitable Earths (APACHE, Sozzetti et al. 2013) and EXOplanetary systems Robotic APT2 Photometry (EXORAP). In the following, we discuss the GLS analysis of this data together with additional photometry from the Wide Angle Search for Planets (WASP, Butters et al. 2010) and Hipparcos (Perryman et al. 1997) catalogs in order to identify at least the rotational periodicity at 16.3 days.

The APACHE photometric survey monitored GJ 3942 with a 40 cm telescope located in the Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA, +45.7895 N, +7.478 E, 1650 m.a.s.l.), between August 30 and September 29, 2013 (12 nights of season 1 of HARPS-N observations), and May 7 to 15, 2014 (4 nights of season 2). The dataset covers a time span of 258 days, and is composed of 418 useful measurements. The observations were carried out using a Johnson V filter following the standard strategy used by APACHE, consisting of three consecutive exposures repeated at intervals of ~20 to 25 min while the target is >35° above the horizon. The images were reduced with the standard pipeline TEEPEE by the APACHE team (see Giacobbe et al. 2012). The light curve, including 141 data points, was obtained by using an aperture of 6.5 pixels (9.75 arcsec) and a set of five comparison stars (UCAC4 716-054601, UCAC4 716-054610, UCAC4 715-054335, UCAC4 715-054346, and UCAC4 715-054359). To analyze the data with the GLS, we calculated nightly averages (see panel a of Fig. 6) with a dispersion of 5.9 mmag and average errors of 4.2 mmag. The GLS periodogram of these 16 data points is shown in panel (a) of Fig. 7. A best fit is found for a period of 18.1 days and a semi-amplitude of 7.0 mmag with a significance quite close to the 1% threshold.

EXORAP is carried out at INAF-Catania Astrophysical Observatory with a 80 cm f/8 Ritchey-Chretien robotic telescope (Automated Photoelectric Telescope, APT2) located at Serra la Nave (+14.973°E, +37.692°N, 1725 m.a.s.l.) on Mt. Etna. BVRI photometry of the star was collected over 85 nights between May 5, 2014, and April 30, 2017. The data cover seasons 2 to 4 of the HARPS-N observations. To obtain differential photometry, we started with an ensemble of ~6 stars, the nearest to GJ 3942 having similar brightness. Then we checked the variability of each of them by building their differential light curves using the rest of the sample as reference. That way we selected the four least variable stars of the sample. The rms of the ensemble stars is 10, 13, 19, 25 mmag in B, V, R, and I, respectively. We obtained 100/89/87/86 data points for the B/V/R/I filters (see panels b, c, d, and e of Fig. 6), respectively, and no data points were rejected by a 5σ clipping. We calculate average dispersions of 15.1/12.7/17.3/20.1 mmag and average photometric uncertainties (sky + poisson) of 1.1/0.9/1.2/1.1 mmag for the four filters. For every dataset, we correct for a long-period periodicity (>250 days). The GLS periodograms of the residuals are shown in Fig. 7 for filters B (panel b), V (panel c), R (panel d), and I (panel e). The 16.3-day signal is seen in every dataset, but it is at the 1% FAP level in the B filter and very close to it in V and R. Best-fit periodicities are therefore 16.34/17.04/16.30 days with semi-amplitudes of 6.5/6.6/7.8 mmag in the BVR filters. Those values reproduce very well our previous estimations of the rotational period. The fits reduce the respective overall variation by 11.5/12.7/16.8%. We do not show the WFs of all the photometric time series data since they do not influence their respective periodograms in the important period range.

We used WASP2 optical photometry containing 7447 data points. The target was observed in 120 nights from April 12, 2006, to April 30, 2008. We show the data after 1σ clipping of the RVs and RV uncertainties and after averaging the values of each night in panel (f) of Fig. 6. It shows a mean magnitude of V = 10.234 ± 0.044 mag and a mean error of 13.2 mmag. No significant periodicities are found in the GLS analysis. A last dataset of 16 nightly averages of 116 measurements from the Hipparcos catalog are shown in panel (g) of Fig. 6. They were observed from January 3, 1990, to March 8, 1993, in the HP filter, which is similar to Johnson’s V. We calculate an average magnitude of 10.242 ± 0.021 mag and an average uncertainty of 12.2 mmag. The GLS analysis of the data does not show any significant signals.

The photometric data of the APACHE and EXORAP programs confirm the 16-day rotational period of GJ 3942. All smaller signals shorter than <100 days are related to this main periodicity and are corrected for after pre-whitening the rotational periodicity. For the EXORAP data, which cover the same seasons as our spectroscopic HARPS-N observations, we can identify a minimum in brightness by season 3 (mid-2015 or JD ~ 2 457 200 days).

thumbnail Fig. 7

GLS periodograms of the different photometric datasets including the APACHE V filter (panel a) and EXORAP filters B (panel b), V (panel c), R (panel d), and I (panel e). The yellow bands indicate the periods at 16.3, 10.4, and 6.9 days. The blue dashed horizontal line indicates the 1% FAP level.

Open with DEXTER

3. Detailed analysis

Our analysis of the available RV time series data of GJ 3942 suggests the presence of a planetary companion with an orbital period of 6.9 days and of another signal at a period of 16.3 days that we attribute to the stellar rotation period as traced by magnetic regions on the stellar surface. In the following, we compare different methods in order to find arguments for those findings and to further investigate whether the remaining signals could be suggestive of the presence of a second planet or are instead related to magnetic activity arising from the evolution and decay of active regions and/or WF aliasing. Here we use the classical iterative method using a modified GLS code, the likelihood-ratio periodograms including a moving average term, and the Gaussian-process regression.

thumbnail Fig. 8

GLS periodograms of the RV data of GJ 3942. The contents of the panels are a) the periodogram of the RV data; b) residuals after correcting for a 16.3-day signal; c) residuals after subtracting simultaneously 16.3- and 6.9-day signals; and d) residuals after simultaneously correcting for signals at 16.3, 6.9, and 10.4 days. For each panel, we give the analytical FAP of the respective model (in red if FAP < 0.1%). Yellow bands indicate the periods at 16.3, 10.4, and 6.9 days. The green lines show the first and second harmonics of the 16.3-day signal, namely 8.2 and 5.4 days. Blue horizontal lines show the 1% (dashed) and 0.1% (solid) FAP levels (from bootstrapping) of the respective dataset. Additionally, we show in the upper parts of each periodogram the WF in gray after being flipped, scaled, mirrored, and period-shifted to the main peak of the respective dataset.

Open with DEXTER

3.1. M-GLS analysis

To further analyze the data, we eliminated the best sinusoidal fit to the most significant signal in the RV time series. We then recurrently continued this procedure with the residuals in an iterative way following a procedure called pre-whitening (e.g., Hatzes 2013). Instead of correcting for each signal one after the other, we used a generalization of the GLS algorithm that considers a simultaneous multi-frequency fit, i.e., a GLS with more than one dimension, called multi-dimensional GLS (M-GLS). In this way, the method employs an approach that mitigates cross-talk between signals having similar frequencies and amplitudes that could affect the classical pre-whitening technique.

thumbnail Fig. 9

Likelihood-ratio periodograms of the RV data of GJ 3942. The contents of the panels are a) the periodogram of the RV data; b) residuals after correcting for a 16.3-day signal; c) residuals after subtracting simultaneously 16.3- and 6.9-day signals; and d) residuals after simultaneously correcting for signals at 16.3, 6.9, and 5.4 days. The rest of the symbols are the same as in Fig. 8.

Open with DEXTER

We show in panel (a) of Fig. 8 the GLS periodogram of the RV data scaled to its main peak. Only the period range from 4.4 to 50 days (0.02 to 0.22 days-1) is plotted because of the relevance to this study and the lack of additional long-period signals. At the top of the panel the WF is shown in gray. It is period-shifted and mirrored to the main peak of the periodogram, and flipped and scaled for display purposes. This way, it shows the influence of the combination of the time-sampling and the selected signal on the periodogram. Most importantly in our case, it highlights yearly aliases seen on both sides of a peak. The sinusoidal fit parameters of the 16.3-day signal are shown in Table 3 under Nsignal = 1, and the GLS periodogram of the residuals, which have their rms reduced by 22%, in panel (b). After pre-whitening the strongest signal, other features at short frequencies decrease in power, but the 6.9-day peak does not vary. In addition, the two peaks at 10.1 and 10.4 days increase their power level and reach a FAP of 1%. There are also significant changes at around 5.4 days, and we hypothesize that this peak rising up to the 1% FAP level can be related to the second harmonic of the 16.3-day signal, arising from the deviation of the rotational modulation effect from a purely sinusoidal shape. Panel (c) in Fig. 8 shows the periodogram after simultaneously correcting for the 16.3-day signal (see exact parameters in Table 3 under Nsignal = 2) and the next most prominent peak at 6.9 days, which we attribute to GJ 3942 b. The rms of the RV data is reduced by a further 12%. This procedure cleans the period at 6.9 days efficiently, but at the same time a yearly alias of the 16.3-day signal gains slightly in power. The 5.4-day signal becomes slightly less significant and the 10.4-day signal, now the most prominent peak, is close to reaching the 0.1% FAP level. The same procedure is used to remove this third additional signal, reducing the RV rms by another 7%. The fit uses the parameters shown in Table 3 with Nsignal = 3. The GLS periodogram of the residuals in panel (d) indicates that no statistically significant signal is left, although the strongest peak appears at 5.4 days. We note that the subtraction of the 10.4-day signal also removes the signal at 10.1 days as they are yearly aliases of one another: (10.09-1–10.38-1)-1 ≈ 365 days. The classical pre-whitening approach using the M-GLS algorithm therefore favors three periodic signals of 16.3, 6.9, and 10.4 days.

Table 3

Best-fitting parameters of the RV data of GJ 3942 and its residuals during the pre-whitening procedure of the three methods applied.

3.2. Likelihood ratio periodograms

To further analyze the data independently, we use likelihood ratio periodograms on our RV dataset. This includes considering Keplerian fits of planetary orbits and periodic stellar signals as well as a first-order moving average (MA) component with exponential smoothing to account for the remaining correlated signals as red-noise (see, e.g., Feroz & Hobson 2014; Tuomi et al. 2013). For the details of the likelihood function, which is calculated by the best fits to the time series data, we refer to Anglada-Escudé et al. (2016). The additional term of the first-order MA accounts for the correlated noise considering the RV residual ϵ of the previous measurements at ti − 1. The MA at each time ti is then given by (3)The coefficient φ measures the correlation of RVs at ti and ti − 1 and τ is a characteristic time decay measuring the impact of that correlation to the RVs over time.

In Fig. 9 we show equivalent panels to those in Fig. 8 using the likelihood ratio periodograms. The plots are again scaled to the main peak of the periodogram in panel (a). In general, thanks to the more sophisticated treatment of data correlations, the background noise level is diminished when using this technique. We show the different parameters used for the fits in Table 3 and the additional parameters used by this technique in Table 4. Panel (a) shows two strong peaks, as before, at approximately 16.3 and 6.9 days, weaker double peaks at around 10 days, and a moderately strong periodic signal at 8.2 days, which is interpreted as the first harmonic of the 16.3-day period. As can be seen in panel (b), removing the 16.3-day signal diminishes the significance of its first harmonic and increases the power of the 6.9-day signal, as do again the signals at approximately 5 days. The RV rms is reduced in this step by 33%. Removing both the 16.3- and 6.9-day signals reduces the rms by another 14%, and the 5.4-day peak, which has some structure, becomes the most significant one, as shown in panel (c). Subtracting the 5.4-day signal leaves some power at a similar period of 5.2 days, as can be seen in panel (d). Also, we have some excess power left at approximately 8, 10, 16, and 20 days. The reduction of the rms by this last step is on the order of 8%.

The MA method seems to favor a one-planet model with a slightly lower FAP when compared to the M-GLS method and small orbital and effective eccentricities, respectively. Together with the red-noise treatment, the method is thereby able to reduce the rms of the residuals significantly in comparison with the M-GLS analysis. The signal at 5.4 days appears strong after correcting for the rotational periodicity in (b) indicating again that it might be the second harmonic of the 16.3-day signal. A further remaining peak at 5.2 days in (d) could then represent a first harmonic of the signal at 10.4 days that we find with the M-GLS method.

Table 4

Additional best-fitting parameters of the RV data of GJ 3942 and its residuals using the MA method.

3.3. Gaussian process regression

thumbnail Fig. 10

Distributions of the total MCMC samples for the jump parameters of the GP one-planet model as a function of the natural logarithm of the likelihood function (yellow dots). The samples corresponding to λ> 100 days are shown as gray dots. The vertical lines represent the median values (solid) and the 16th and 84th percentiles (dashed) of the latter distribution, which are listed in Table 5.

Open with DEXTER

Following our multi-technique approach, we modeled the RVs of GJ 3942 by using Gaussian process (GP) regression through a Markov chain Monte Carlo (MCMC) algorithm. This method is becoming increasingly popular as an efficient and physically robust way to mitigate signals in RV data that can be ascribed to stellar activity while retrieving reliable estimates of the planetary orbital parameters. For details of how GPs have been used for the analysis of RV datasets, we refer to, e.g., Haywood et al. (2014), Grunblatt et al. (2015), Rajpaul et al. (2015), Affer et al. (2016), Faria et al. (2016), López-Morales et al. (2016), Dumusque et al. 2017, Damasso & Del Sordo (2017), Cloutier et al. (2017), Angus et al. (2017). In this work we have adopted the widely used quasi-periodic kernel described by the following covariance matrix: (4)where t and t indicate two different epochs. This functional form is particularly suitable to modeling periodic short-term stellar signals, modulated by the stellar rotation period, by allowing for an exponential decay of the correlations between different epochs. The hyper-parameters appearing in Eq. (4) are linked to some of the physical phenomena underlying the stellar noise: h represents the amplitude of the correlations, θ generally represents the rotation period of the star, w is the length scale of the periodic component linked to the size evolution of the active regions, and λ is the correlation decay timescale that we assume to be related to the active regions lifetime. Moreover, σRV(t) is the RV internal error at time t, σjit is the additional uncorrelated “jitter” term that we add in quadrature to the internal errors to take into account instrumental effects and other noise sources not included either in σRV(t) or in our stellar activity framework, and δt,t is the Dirac delta function.

We refer to Damasso & Del Sordo (2017) for the log-likelihood function to be maximized by the MCMC procedure and the general form for the models that we tested in this work. The selection of the model, multiplicity, and the priors of the GP hyper-parameters are based on the outcome of the M-GLS and MA analyses of the previous sections.

3.3.1. One-planet model

We list in Table 5 the uniform prior probability distributions for the parameters of the one-planet model. Regarding the GP hyper-parameters, w was constrained between 0 and 1 from the discussion in López-Morales et al. (2016), θ was sampled in a range around the stellar rotation period of 16.3 days, and an upper limit for λ was adopted with a value slightly higher than the dataset timespan. The planetary orbital period P1 was sampled in a small interval around 6.9 days. The parameter space was explored using 150 random walkers. We randomly initialized the positions for the λ GP hyper-parameter around 100 days using a Gaussian distribution (σ = 50 days) to start the exploration of the parameter space within a sufficiently wide range for the correlation decay timescale. The progress of the MCMC fitting procedure was monitored by evaluating the Gelman-Rubin convergence parameters as defined by Ford (2006). The final distributions of the free parameter samples as a function of ln(likelihoodprior) are shown in Fig. 10 by the yellow dots. These were obtained by applying a first burn-in of 3000 steps and then discarding additional samples up to the first MCMC step at which all chains have had at least one value of ln(likelihoodprior) lower than the median of the ln(likelihoodprior) dataset, following the prescription of Eastman et al. (2013, and references therein). We note that the posterior distribution of the λ timescale appears to be bi-modal. The maximum a posteriori probability estimate λ = 17.1 days, which is the expected peak at about the stellar rotation period. Of astrophysical relevance is the second local maximum at ~250 days (Strassmeier 2009; Davenport 2015). To derive the final best-fit parameter values, we then consider only samples for which λ> 100 days (the period where both solutions meet), resulting in the posterior distributions shown in Fig. 10 by the gray dots. The uncorrelated jitter appears to be bimodal, whereas for our 250 day λ the solution occurs at ~2.7 m s-1, which is more than double the mean RV error. The high value for the additional jitter could be indicative of the signal of a second planetary in the dataset. Table 5 shows the corresponding best-fitting values and uncertainties for each jump parameter, calculated as the median of the marginal posterior distributions and the 16% and 84% quantiles. In Table 3 we show the characteristics of the fit to compare them to the other methods. We note that the θ parameter from the GP term is regarded as a period of the first signal, but given in parentheses.

The values of the one-planet model reproduce the results of the other methods very well. The value for λ confirms the described changes of magnetic activity phenomena on a seasonal timescale. Figure 11 shows the stellar contribution to the RV after removing the best-fit planetary solution. We note a slight increase in the scatter over the four seasons, corresponding to the RV values. The residuals of our global model show an rms of 2.36 m s-1. The Bayesian strength of the models tested in this work are evaluated by using the Bayesian information criteria (BIC), and the results are interpreted by adopting the empirical scale presented in Raftery (1995). For the one-planet model we obtain BIC = 861.3.

thumbnail Fig. 11

Radial velocity residuals time series (black dots) after subtracting our best-fit orbital solution for GJ 3942 b. The blue line with gray shaded 1σ regions represents our best-fit GP quasi-periodic model for the correlated stellar noise. The top plot shows the complete dataset, while the lower plot shows a blow-up of the third epoch for easier visualization of the agreement between the model and the data.

Open with DEXTER

Table 5

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters used in the one-planet circular model.

3.3.2. Modeling activity index time series

The same approach was used to analyze the time series of the S and Hα indices. We performed a GP regression without including an additional jitter term in the kernel, and loosely conditioning the model a priori. The hyper-parameter θ was constrained within a small range around a 16.3-day period, while we left λ free. The positions of the random walkers were initialized at around 550 days. The uniform priors are listed in Table 6. The 50 independent chains reached convergence according to the criteria defined by Ford (2006). We show in Fig. 12 the posterior distributions obtained after a burn-in as explained in the previous section. We see that the stellar rotation period θ appears to be slightly bi-modal but still tightly constrained around 16.3 days, thus reproducing the solution found fitting a simple sinusoid. The best-fit solution has very similar properties to that found for the stellar contribution in the RVs alone, but with a larger λ of 634 days. Table 6 lists the corresponding best-fitting values and uncertainties for each jump parameter.

Table 6

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters of the quasi-periodic kernel used to model the S index dataset.

thumbnail Fig. 12

Posterior distributions of the GP hyper-parameters for the quasi-periodic kernel applied to the S index time series. Dashed vertical lines mark the median values and the 16th and 84th percentiles of these distributions, which are listed in Table 6.

Open with DEXTER

Applying the quasi-periodic kernel to the Hα index dataset does not produce the same results as for the S index. By adopting the same priors, the chains did not reach convergence with θ = 16.262 ± 0.007 days and λ peaking at the upper edge of the prior interval. This indicates very stable fit parameters over the time of observations. The result is in agreement with the lowest FAP of this index in Fig. 3 for the rotational signal, which indicates that it is best described by a simple sinusoid. It also agrees with the smallest variation in the fit parameters of season 3 in Fig. 4.

3.3.3. Two-planet model

As discussed before, the analysis of the RV residuals suggests the existence of an additional signal. To explore the significance of this possibility, we tested a model with two circular Keplerian signals by using the same uniform priors of the one-planet model, and a third periodicity between 0.1 and 25 days (see Table 7). This time we used 100 independent MCMC chains, initially spread over a large fraction of the parameter space, and we stopped the run after 200 000 steps without reaching the formal convergence. The results show that the posterior distribution for P2 is characterized by two local maxima at approximately 10.1 and 20.4 days, with the latter period being more significant. The solution points towards an evolutionary timescale λ that is similar to the one-planet model, and a slightly smaller semi-amplitude for the inner planet. The residuals have an rms of 1.70 m s-1, well below the expected jitter of 2.3 m s-1. In the residuals the 10.1-day periodicity is still visible. The BIC value for this model is 856.7. By using the approximation BIC2 − BIC1 ≃ 2lnB12, we obtain B12 ~ 0.1 for the Bayes factor, indicative of mildly positive evidence for the two-planet model. Models with free eccentricity did not converge.

Table 7

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters used in the two-planet circular model.

4. Discussion

Two main periodic signals in our data are confirmed by all the methods that we employ, one that is related to stellar magnetic activity and one that is best described as induced by the Keplerian motion of an orbiting planet. The differences between the various analysis methods appear in the strength and interpretation of subsequent signals, most precisely of a third RV periodicity that might be an additional planet in the system (see Table 3).

The M-GLS fit proposes such a second planetary candidate with an orbital period of 10.4 days. We model rotational and planetary modulations as simple sinusoids with constant parameters, although they might have evolving characteristics and eccentric orbits, respectively. The method also does not treat the correlated or uncorrelated noise components, which make this simple correction that uses the pre-whitening approach more unreliable. The likelihood-ratio MA analysis instead includes a noise contribution into the fits. It favors an interpretation of harmonics of the main signal at 16.3 days (P) instead of an additional independent signal. This is seen in panel (a) of Fig. 8 for the 6.9-day and P/2 (8.15 days) signals and in panel (c) for the 10.4-day and P/3 (5.4 days) signals. Thus, the not significant 5.4-day peak selected as a third periodicity by this method would be a harmonic of the rotation period and the presence of such harmonics would be a consequence of the varying and non-sinusoidal shape of the 16.3-day modulation in the RV time series data. Also, we suspect some correlation between the time sampling and the periodic signals found since we can reproduce the 5.4 day peak by simple models simulating a dataset with mildly eccentric 16.3, 10.4, and 6.9-day Keplerian signals. The GP regression method adjusts the observed data to a certain model of the stellar system and results in the detection of a 20.4-day period in addition to the strong 16.3- and 6.9-day signals. The suggestive harmonic chain at 20.4, 10.2, 6.8, and 5.1 days, which are close to some of the significant signals that we have found throughout our study, deserves further analysis.

In the case of an active star, the presence of differential rotation, i.e., latitudinal differences in the rotation velocity, can manifest itself through the appearance of additional signals in the periodogram. The analysis of Kepler data (Balona et al. 2016) shows that differential rotation is common in all stars, but especially so in slow and moderate rotators such as GJ 3942. For example, Reinhold & Gizon (2015) used the α parameter (α = 1 − Pmin/Pmax) to quantify the level of differential rotation through the empirical analysis of 12 000 stars with photometry from the Kepler mission. We hypothesize that the 20.4-day signal revealed by the GP method arises from active regions at high-latitude locations. Indeed, the α parameter of this period compared with 16.3 days is 0.201, in agreement with Fig. 9 of Reinhold & Gizon (2015). There is, however, a strong indication that the 20.4-day signal and its harmonics are not associated with stellar activity. The reason is the clean periodogram of the activity indicators S and Hα, where only the 16.3-day signal and, possibly, its much weaker first and second harmonic signals are seen. There are no traces of signals at 20.4 days or any of its harmonics. Of course, using this argument to rule out the 20.4-day harmonic chain as arising from stellar activity hinges on the fact that the chromospheric emissions traced by these indicators are good proxies for photospheric active regions that produce the RV shifts. This seems to be a reasonable assumption as such a connection is observed in the Sun and we note that we see only the 16.3-day periodicity in the photometric datasets.

thumbnail Fig. 13

Residuals and errors of the RVs in black dots folded to a period of 20.5 (top panel) and 10.4 days (bottom panel) with semi-amplitudes of 1.88 and 2.52 m s-1, respectively. We used the data after pre-whitening the 16.3-day signal (top) and the 16.3- and 6.9-day signal (bottom) using the MA method. The red dots correspond to averages in 0.05 phase intervals. We show the best sinusoidal fits in blue.

Open with DEXTER

Another scenario could be that the harmonic chain is produced by a very eccentric planet. For illustration, we show in the top panel of Fig. 13 the RV residuals of our MA method after correcting for the periodic signal at 16.3 days and phase folded at a period of 20.52 days, which is the best-fit period in an interval from 20 to 21 days. The resulting sinusoidal modulation has a semi-amplitude of only 1.88 m s-1, which is well below the expected uncorrelated RV jitter of 2.3 m s-1, and does not show any signs of high eccentricity. We also carried out a simulation analysis using the observed time series dates and injecting signals at 16.3 and 20.4 days and considering Gaussian errors. After performing 10 000 trials by varying eccentricities for the rotational and planetary signals from 0 to 0.3 and from 0.3 to 0.9, we did not find any configuration that could reasonably reproduce the set of signal harmonics. From this, we deem the 20.4-day period and its harmonic chain as quite unlikely to arise from a highly eccentric planet, or actually by more than one since at these close periods the system would be dynamically unstable.

thumbnail Fig. 14

Evolution of periodogram power of the signal at around 10.4 days with the number of observations. We use the residuals after the subtraction of the 16.3- and 6.9-day signals of the best fit following our MA method in red and in the original data in black. The solid and dashed lines correspond to the 0.1 and 1% FAP levels, respectively.

Open with DEXTER

Table 8

Summary of the characteristics of the RV periodicities of GJ 3942 following the MA method.

A further interesting hypothesis is that the signals are induced by a true planetary system of two or more planets in near resonant orbits. Orbital resonances are found to be quite common in exoplanet systems, including the TRAPPIST-1 system, with up to seven planets near mean-motion resonances (Gillon et al. 2017; Luger et al. 2017). The strengths of the various signals of the GJ 3942 RVs are still too weak to make this claim, but it is worth studying the most prominent of the periods resulting from the frequentist analysis at 10.4 days. Its period ratio with GJ 3942 b at 6.9 days is 3:2. This could be a strong sign of a 3:2 resonance which is, in contrast to the 3:1 resonance the 20.4-day signal would represent, a quite common configuration in multi-planetary systems as found by the Kepler mission (Lissauer et al. 2011; Wang & Ji 2014). We show the chronological evolution of the periodogram power of the 10.4-day signal in Fig.14. We calculate the evolution both in the original data and also after subtracting the 16.3- and 6.9-day periods. As expected from a nonspurious signal, the significance improves with the treatment of the RV data. The signal came close to reaching the 0.1% FAP level, but the last ten observations have somewhat degraded the significance. In the bottom panel of Fig. 13 we show the RV residuals of our MA method, after correcting for the periodic signals at 16.3 and 6.9 days and phase-folded with a period of 10.38 days. The resulting fit has a semi-amplitude of 2.52 m s-1. To address the dynamical stability of the putative system, we use the definitions given by Giuppone et al. (2013) and their Eq. (3). As a result, we find that a two-planet system with orbital periods of 6.9 and 10.4 days and masses of 7.2 and 6.3 M can be dynamically stable, albeit close to the limit and only for nearly circular orbits (eccentricity < 0.2).

After considering the various scenarios, we can only confidently claim the existence of one planet around GJ 3942, which we dub GJ 3942 b and which is likely to be a super-Earth with a minimum mass of 7.1 M in a mildly eccentric orbit with a period of 6.9 days. Our best Keplerian fit of the periodicity with the MA method is shown on the bottom panel of Fig. 15, the detailed characteristics are given in Table 8. To be able to estimate uncertainties of the fitted parameters of the Keplerian orbit (Eq. (2)), we use the IDL routine RVFIT (Iglesias-Marzoa et al. 2015) on the respective RV time series data with the limits for the values as shown in the table. A second planet in the system cannot be confirmed from our analysis of the available data, and we thus deem it to be a candidate, with its possible properties also listed in Table 8. The rotational modulation signal is the strongest in the RV time series data (and also in the activity indices and photometry) and is shown in Fig. 15 (top) and Table 8 if modeled using the same parameters as Keplerian motion.

thumbnail Fig. 15

Best simultaneous Keplerian fits (blue curve) following our MA method for periods of 16.28 (top) and 6.91 days (bottom). The fits have semi-amplitudes of 5.3 and 3.3 m s-1 and eccentricities of 0.20 and 0.12, respectively. The data are shown without the contribution from the MA term.

Open with DEXTER

5. Conclusions

We investigated 145 spectroscopic observations of GJ 3942 obtained over 3.3 yr with HARPS-N at the TNG in La Palma, Spain, and additional photometry from various sources, in particular the APACHE and EXORAP programs. We used RVs from the TERRA pipeline and activity indices calculated from Hα and Ca ii emission lines originating in the stellar chromosphere. We investigated further the instrumental drifts of the observations and proposed adding quadratically around 1 m s-1 to the RV uncertainty if the second HARPS-N fiber is not used for simultaneous calibration (e.g., due to long exposure times).

Three different approaches were used to search for planets. Significant in all of them is the one-planet model for which the moving average method delivers the lowest false alarm probability. Therefore, and since it includes a treatment of the red noise – in contrast to the multi-dimensional generalized Lomb-Scargle periodogram – the results of this technique are the most reliable. The detailed Gaussian process regression, although probably best able to reproduce the changing characteristics of the magnetic periodicity, cannot resolve the GJ 3942 system properly since it adds the not significant but existent 10.4-day signal into the noise or Gaussian process term. The one-planet model consists of two significant signals confirmed by all our considerations.

A dominant periodic signal at 16.28 days is clearly visible in the photometric data, the activity indices, and the RV measurements. The RV data and the chromospheric indices show a phase shift of 120 deg and we observe changes in the characteristics of the magnetic activity signal on rotational to seasonal timescales. We conclude that this period is closely related to the rotation of GJ 3942 and the evolution of magnetically active regions on the stellar surface. Because of the complex nature, those signals cannot be described by single sinusoids. Together with possible resonances between the time-sampling and present periodicities, this introduces power at 8.15 and 5.4 days, which are the first and second harmonics, respectively. It is interesting to note a rather sharp variation in the magnetic activity of the star at the end of the third season (August 2015), which also seems related to a brightness minimum in the photometric EXORAP data. These variations can also be connected with the decreasing strength of the activity indices and the increasing RV scatter from seasons 1 to 4. The pattern changes also affect the detectability of further periodicities.

The second strongest signal in the periodogram at 6.91 days has a high significance level that is not diminished by the modeling of magnetic activity. This periodic signal is not present in any activity index nor the photometric data and has a very stable period value. We therefore propose it to be induced by the orbit of a planet, GJ 3942 b. From our best fit we find a mass of 7.14 ± 0.59 M and an orbital semi-major axis of a = 0.0608 ± 0.0068 AU (see Table 8). Using zero albedo, we calculate roughly an equilibrium temperature of 590 K using .

Furthermore, we find still inconclusive evidence of at least one more planetary signals in our data, at 10.38 days. Its significance rose steadily over the observational time span and it is prominent in our analysis, but still lacks significance. The periodogram signal shows a visible harmonic at 5.2 days and a yearly alias at 10.1 days. This second planet, GJ 3942 c, would orbit with a 3:2 period ratio with GJ 3942 b and would have an equilibrium temperature of 515 K. The possibility of additional planets in the system, likely in a near-resonant chain, is tantalizing but not confirmed by our analysis, and thus remains open.

It is also interesting to assess the possibility of detecting transits from GJ 3942 b. We estimate the geometric transit probability at 4.7% and a maximum duration of 2.5 hours. From the stellar rotational period of 16.28 ± 0.22 days and a radius of 0.61 ± 0.06 R we calculate a rotational velocity of v = 1.89 ± 0.19 km s-1. If we compare this with the vsini measured from spectral features by Maldonado et al. (2017), 1.67 ± 0.30 km s-1, we find it to be compatible with an inclination of the rotational axis of 90 deg. Adopting the mass-radius relationship of exoplanets described in Fig. 8 of Dumusque et al. (2014), and assuming a 1.8 R for a rocky configuration and 3 R for a low-density planet, we calculate possible transit depths of 0.8 to 2.2 mmag. This is challenging to detect with ground-based telescopes, but could be an excellent target for the upcoming CHEOPS mission. The HADES observations on the other hand, are still ongoing in order to be able to find more hints of the existence of further planets around GJ 3942, and to fulfill the prime goal of the program, which is the discovery of a statistically significant sample of M-dwarf planets to provide stronger constraints on their occurrence rates, system architectures, and formation mechanisms.


Acknowledgments

We are grateful to Guillem Anglada-Escudé for the HARPS-N TERRA code and the Likelihood-ratio diagrams. M.P., I.R., J.C.M., A.R., E.H, and M.L. acknowledge support from the Spanish Ministry of Economy and Competitiveness (MINECO) and the Fondo Europeo de Desarrollo Regional (FEDER) through grant ESP2016-80435-C2-1-R, as well as the support of the Generalitat de Catalunya/CERCA program. L.A. and G.M. acknowledge support from the Ariel ASI-INAF agreement No. 2015-038-R.0. A.S.M., J.I.G.H., R.R.L., and B.T.P. acknowledge financial support from the Spanish Ministry project MINECO AYA2014-56359-P. J.I.G.H. also acknowledges financial support from the Spanish MINECO under the 2013 Ramón y Cajal program MINECO RYC-2013-14875. GAPS acknowledges support from INAF through the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement No. 313014 (ETAEARTH). The HARPS-N Project is a collaboration between the Astronomical Observatory of the Geneva University (lead), the CfA in Cambridge, the Universities of St. Andrews and Edinburgh, the Queen’s University of Belfast, and the TNG-INAF Observatory. This paper makes use of data from the first public release of the WASP data as provided by the WASP consortium and services at the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has made further use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

All Tables

Table 1

Stellar parameters of GJ 3942.

Table 2

Different methods used to estimate instrumental drifts and the resulting rms of the data.

Table 3

Best-fitting parameters of the RV data of GJ 3942 and its residuals during the pre-whitening procedure of the three methods applied.

Table 4

Additional best-fitting parameters of the RV data of GJ 3942 and its residuals using the MA method.

Table 5

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters used in the one-planet circular model.

Table 6

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters of the quasi-periodic kernel used to model the S index dataset.

Table 7

Uniform prior probability distributions and best-fitting estimates for the hyper-parameters used in the two-planet circular model.

Table 8

Summary of the characteristics of the RV periodicities of GJ 3942 following the MA method.

All Figures

thumbnail Fig. 1

All available 8199 instrumental drifts distributed over 356 nights. The solid line indicates the mean value and the dashed lines the 1σ interval.

Open with DEXTER
In the text
thumbnail Fig. 2

Time series data of (from top to bottom) radial velocities (in m s-1), S index, and Hα index (in 10-2) are shown in panels a, c, and e, respectively, as observed (black dots) and their linear trend is indicated by the blue lines. The mean errors of each dataset are illustrated by the red line on the lower left side of each panel. Panels b, d, and f show the de-trended data (i.e., linear fit subtracted) folded to the most prominent periodicity at 16.3 days and with zero phase corresponding to the first observation (black dots) and a simple sinusoidal fit (blue curve).

Open with DEXTER
In the text
thumbnail Fig. 3

Generalized Lomb-Scargle periodograms of the de-trended time series data and of the residuals after correcting for the 16.3-day signal with a simple sinusoidal fit. In panels a and b this is shown for the radial velocities, in panels c and d for the S index, and in panels e and f for the Hα index. Additionally, we show in panel g the window function of the time series. Yellow lines indicate the important periodic signals for GJ 3942, namely 16.3, 10.4, and 6.9 days, and we also show additional periodicities mentioned in the text. The dashed blue horizontal line shows the 1% FAP level of the respective dataset.

Open with DEXTER
In the text
thumbnail Fig. 4

Illustration of changing characteristics of signals in RV, RV residual, and activity index data of GJ 3942. It is grouped into seasons 1 & 2 (S12, green dots), season 3 (S3, red dots), and season 4 (S4, black squares). The best purely sinusoidal fit for each dataset is shown. Panel a shows the RV data folded in phase with a 16.29-day period; panel b the residuals after correcting for the best sinusoidal 16.3-day period fit and phase folded with 6.91 days; panel c the S index data folded in phase with a 16.26-day period; and panel d the Hα index data phase folded with 16.25 days.

Open with DEXTER
In the text
thumbnail Fig. 5

FAP development for the 16.3-day (black) and the 6.9-day (red) signals of the RV data using an increasing number of data points. The dashed horizontal line shows the analytical 0.1% FAP level, and the two vertical dash-dotted lines divide the data into the different seasons.

Open with DEXTER
In the text
thumbnail Fig. 6

Photometry of GJ 3942 including normalized data from APACHE (panel a), magnitude differences (GJ 3942 and calibration stars) of EXORAP filters B (panel b), V (panel c), R (panel d), and I (panel e), and visual magnitudes from the WASP (panel f) and Hipparcos catalogs (panel g). The red lines on the lower left show the rms of the nightly averages (panels a, f, g) and the RV rms of the calibration stars (b to e) as proxies for the uncertainties of the measurements. We note that APACHE was used during seasons 1 and 2 of the HARPS-N observations, EXORAP from season 2 to 4 and beyond, and WASP and Hipparcos data are from around 2007 and 1991, respectively.

Open with DEXTER
In the text
thumbnail Fig. 7

GLS periodograms of the different photometric datasets including the APACHE V filter (panel a) and EXORAP filters B (panel b), V (panel c), R (panel d), and I (panel e). The yellow bands indicate the periods at 16.3, 10.4, and 6.9 days. The blue dashed horizontal line indicates the 1% FAP level.

Open with DEXTER
In the text
thumbnail Fig. 8

GLS periodograms of the RV data of GJ 3942. The contents of the panels are a) the periodogram of the RV data; b) residuals after correcting for a 16.3-day signal; c) residuals after subtracting simultaneously 16.3- and 6.9-day signals; and d) residuals after simultaneously correcting for signals at 16.3, 6.9, and 10.4 days. For each panel, we give the analytical FAP of the respective model (in red if FAP < 0.1%). Yellow bands indicate the periods at 16.3, 10.4, and 6.9 days. The green lines show the first and second harmonics of the 16.3-day signal, namely 8.2 and 5.4 days. Blue horizontal lines show the 1% (dashed) and 0.1% (solid) FAP levels (from bootstrapping) of the respective dataset. Additionally, we show in the upper parts of each periodogram the WF in gray after being flipped, scaled, mirrored, and period-shifted to the main peak of the respective dataset.

Open with DEXTER
In the text
thumbnail Fig. 9

Likelihood-ratio periodograms of the RV data of GJ 3942. The contents of the panels are a) the periodogram of the RV data; b) residuals after correcting for a 16.3-day signal; c) residuals after subtracting simultaneously 16.3- and 6.9-day signals; and d) residuals after simultaneously correcting for signals at 16.3, 6.9, and 5.4 days. The rest of the symbols are the same as in Fig. 8.

Open with DEXTER
In the text
thumbnail Fig. 10

Distributions of the total MCMC samples for the jump parameters of the GP one-planet model as a function of the natural logarithm of the likelihood function (yellow dots). The samples corresponding to λ> 100 days are shown as gray dots. The vertical lines represent the median values (solid) and the 16th and 84th percentiles (dashed) of the latter distribution, which are listed in Table 5.

Open with DEXTER
In the text
thumbnail Fig. 11

Radial velocity residuals time series (black dots) after subtracting our best-fit orbital solution for GJ 3942 b. The blue line with gray shaded 1σ regions represents our best-fit GP quasi-periodic model for the correlated stellar noise. The top plot shows the complete dataset, while the lower plot shows a blow-up of the third epoch for easier visualization of the agreement between the model and the data.

Open with DEXTER
In the text
thumbnail Fig. 12

Posterior distributions of the GP hyper-parameters for the quasi-periodic kernel applied to the S index time series. Dashed vertical lines mark the median values and the 16th and 84th percentiles of these distributions, which are listed in Table 6.

Open with DEXTER
In the text
thumbnail Fig. 13

Residuals and errors of the RVs in black dots folded to a period of 20.5 (top panel) and 10.4 days (bottom panel) with semi-amplitudes of 1.88 and 2.52 m s-1, respectively. We used the data after pre-whitening the 16.3-day signal (top) and the 16.3- and 6.9-day signal (bottom) using the MA method. The red dots correspond to averages in 0.05 phase intervals. We show the best sinusoidal fits in blue.

Open with DEXTER
In the text
thumbnail Fig. 14

Evolution of periodogram power of the signal at around 10.4 days with the number of observations. We use the residuals after the subtraction of the 16.3- and 6.9-day signals of the best fit following our MA method in red and in the original data in black. The solid and dashed lines correspond to the 0.1 and 1% FAP levels, respectively.

Open with DEXTER
In the text
thumbnail Fig. 15

Best simultaneous Keplerian fits (blue curve) following our MA method for periods of 16.28 (top) and 6.91 days (bottom). The fits have semi-amplitudes of 5.3 and 3.3 m s-1 and eccentricities of 0.20 and 0.12, respectively. The data are shown without the contribution from the MA term.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.