© ESO 2019

1 Introduction

Time-series observations with high-resolution spectrographs such as the High Accuracy Radial velocity Planet Searcher of the southern (HARPS; Mayor et al. 2003) and northern (HARPS-N; Cosentino et al. 2012) hemispheres, the Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical echelle Spectrographs (CARMENES; Quirrenbach et al. 2018), or the iodine-cell HIgh Resolution echelle Spectrograph (HIRES; Vogt et al. 1994) are used to detect and confirm planetary companions of stars by the Doppler shifts of their spectra. The measured radial velocities (RVs) show variations over time that are induced by the Keplerian orbit of the planet, or planets, and a contribution of the stellar surface. The latter can include actual motions on the surface such as oscillations or surface granulation with different short-term timescales (Dumusque et al. 2011), but also phenomena such as dark spots and bright plages, which introduce RV variations by reducing the number of photons from either the blue- or the red-shifted side of the rotating host star. Those phenomena are caused by the magnetic field of the star and its detailed structure and variability. The typical lifetime of those surface phenomena follows a parabolic decay law depending on size (Petrovay & van Driel-Gesztelyi 1997) and can lastfor up to several rotation periods (Bradshaw & Hartigan 2014; Scandariato et al. 2017) or even longer (see, e.g., Davenport et al. 2015). In the Sun, we also observe a long-term magnetic cycle of approximately 11 yr¹ that is induced by the migration of the spot patterns across different latitudes and by the variation of the spot number. Such cycles are expected for all stellar types (Lovis et al. 2011; Suárez Mascareño et al. 2016). The stellar contribution to the RVs, therefore, can be described partly as uncorrelated noise, but also possibly correlated with the rotational period of the star, the lifetime of the activity phenomena, or any long-term magnetic cycle.

The disentangling of activity-induced variations and planetary signals is of great interest, especially due to the numerous planet searches around M dwarfs that are currently being conducted. The rotational periods of these cool stars range from some days to months and coincide with the orbital periods of planets in hot or temperate orbits which are, therefore, prime targets for detection and characterization (Kiraga & Stepien 2007; Newton et al. 2016; Díez Alonso et al. 2019). Planets of a given mass and period induce larger RV amplitudes in low-mass stars but these stars are also magnetically more active (e.g., Testa et al. 2015), with longer-living surface phenomena (Giles et al. 2017) and less influence of surface granulation than stars of earlier type (Berdiñas et al. 2016). They are the most numerous stars in the Galaxy and in the immediate solar vicinity (75% of all stars closer than 10 pc are M dwarfs according to the RECONS survey²), giving us the possibility to probe our neighboring planetary population. Most importantly, the surveys around M dwarfs are the best possibility to date to access the domain of rocky planets in order to study their atmospheres and to apply first statistical calculations (e.g., Mayor et al. 2011; Dressing & Charbonneau 2013; Bonfils et al. 2013; Perger et al. 2017a).

Such surveys include the HARPS-N red Dwarf Exoplanet Survey (HADES; Affer et al. 2016), which monitored regularly around 80 M0- to M3-type stars and detected and confirmed seven Neptune- to Earth-like planets (Suárez Mascareño et al. 2017a; Perger et al. 2017b; Pinamonti et al. 2018; Affer et al. 2019). The HADES program is a collaboration between Italian and Spanish institutes, and has also extensively explored the rotational and magnetic behavior of those low-mass stars (Maldonado et al. 2017; Scandariato et al. 2017; Suárez Mascareño et al. 2018; González-Álvarez et al. 2019). Another survey is conducted with the CARMENES instrument, which is monitoring regularly more than 300 M-type dwarfs selected from the input catalog Carmencita (Caballero et al. 2016a). The program has already confirmed a number of planets (Trifonov et al. 2018; Sarkis et al. 2018) and detected up to seven previously unknown planets (Reiners et al. 2018a; Kaminski et al. 2018; Luque et al. 2018; Nagel et al. 2019), including a cold companion of the nearby Barnard’s Star (Ribas et al. 2018). The instrument is an effort of German and Spanish institutes to fully explore for the first time a large and complete sample of nearby M dwarfs of all spectral subtypes. With its wide range of wavelengths, CARMENES is especially capable of disentangling different sources of RV variability. Whereas a planetary orbit shifts the stellar light of all wavelengths equally, the contribution from surface phenomena is connected to physical processes that usually depend on wavelength.

In this work we have combined observations from both the HADES and CARMENES programs in order to search for planetary companions around Gl 49. This is a low-mass star with spectral type M1.5 V located in the Cassiopeia constellation. Its kinematics and large proper motions suggest a membership to the young disk population of the Milky Way (Cortés-Contreras 2016; Maldonado et al. 2017). The star shows a moderate activity level and Sun-like metallicity, and was found to rotate with a period of 18.4 ± 0.7 d by HARPS-N activity indices (Suárez Mascareño et al. 2018) and 19.9 ± 0.4 d by MEarth photometry (Díez Alonso et al. 2019). We provide an overview of its basic properties in Table 1. Gl 49 forms the wide binary system WNO 51 (Washington Double Star catalog; Mason et al. 2001) with a fainter proper-motion companion at 293.1 arcsec to the east, namely Gl 51 (V388 Cas, Karm J01033+623). The companion is a flaring M5.0 V star with a relatively large amplitude of RV variations due to activity (Alonso-Floriano et al. 2015; Tal-Or et al. 2018; Jeffers et al. 2018), but for which reliable membership in young moving groups, and thus an age estimate, has not been reported (e.g., Gagné et al. 2015).

In Sect. 2 we present the spectroscopic observations of Gl 49 and their treatment. In Sect. 3 we analyze the variation and evolution of the activity of the star with additional photometric observations and various time-series data derived from the observed spectra. We study their properties by investigating their periodicities, calculating phase shifts and model the data as correlated noise in a Gaussian process (GP) framework. We describe further in Sect. 4 how we apply different state-of-the-art models to our system in order to ascertain the existence of a planet, and to find the most likely parameters for the proposed system. A detailed discussion about our findings regarding the discovered system and the evolution of the stellar activity over the last six years is given in Sect. 5. We conclude the work in Sect. 6.

2 Spectroscopic observations

We obtained 137 RVs from optical spectra of the HADES program. They were observed over six seasons (S1–S6) between 3 Sep 2012 and 11 Oct 2017 with HARPS-N. The instrument is installed since 2012 at the 3.58 m Telescopio Nazionale Galileo (TNG) located at the Roque de Los Muchachos Observatory in La Palma, Spain, and is connected to its Nasmyth B focus through a front-end unit. It is a fiber-fed, cross-dispersed echelle spectrograph with a spectral resolution of 115 000, covering a wavelength range from 3830 to 6900 Å. In the HADES framework, we observed all targets with fixed integration times of 900 s to obtain data of sufficient signal-to-noise ratio(S/N) and to average over short-term variations. The data were reduced with the data reduction pipeline DRS (Lovis & Pepe 2007). We extracted RVs using the Java-based Template-Enhanced Radial velocity Re-analysis Application code (TERRA; Anglada-Escudé & Butler 2012), which has been shown to deliver RV time-series of lower dispersion than the binary mask technique used by the DRS, at least in the case of early M-type stars (Perger et al. 2017a). As outlined by those authors, we added quadratically a value of 1.0 m s⁻¹ to the RV uncertainty if no correction for RV drifts was applied by observing simultaneously a Th-Ar lamp with the second fiber. The DRS pipeline delivers a value for the absolute RV shift for Gl 49 of −5.777 ± 0.066 km s⁻¹. The statistics of the datasets are shown in Table 2, and the RV values are given in Table A.1 and visualized in the top panels of Fig. 1. The whole HARPS-N set shows small RV uncertainties but a strong variation in RV scatter over the seasons. Assuming a semi-periodic behavior, we would set a lower limit of recurring apparent RV amplitudes at approximately 1500 d, with a minimum at S3–S4. Additionally, S1 shows an exceptionally large RV scatter, pointing toward a rather variable stellar contribution.

Furthermore, we obtained spectroscopic observations with the CARMENES instrument, installed since 2015 at the 3.51 m telescope of the Calar Alto Observatory in Spain. Its wavelength coverage ranges from 5200 to 9600 Å in the optical and up to17 100 Å in the near-infrared, with resolutions of 94 600 and 80 400, respectively (Quirrenbach et al. 2016, 2018). The visual channel is capable of reaching RV precisions of 1–2 m s⁻¹ (e.g., Trifonov et al. 2018). As reduction pipeline, we used the CARMENES Reduction And CALibration software (CARACAL; Caballero et al. 2016b), and to calculate the RVs, the template-based SpEctrum Radial Velocity AnaLyser (SERVAL; Zechmeister et al. 2018). In the framework of the CARMENES Guaranteed Time Observation program, the target was observed 80 times from 7 Jan 2016 to 25 Feb 2018, covering the three seasons S4, S5, and S6 of the HARPS-N observations. The optical echelle spectra were obtained with exposure times from 100 to 800 s in order to reach a S/N of 150. The optical data set shows a low RV scatter. The data of the near-infrared channel of CARMENES, on the other hand, are not used in this study since the RV amplitudes that we are studying for our early M-dwarf target are expected to be smaller than the RV precision.

Gl 49 was also observed with the HIRES instrument, installed since the late 1990s at the Keck I telescope located in Hawaii, USA (Vogt et al. 1994). The instrument observes from 3000 to 10 000 Å with a resolution of up to 85 000 and a precision of down to 1 m s⁻¹. HIRES uses a iodine cell to monitor and correct for the RV drifts introduced by temperature changes. In total, data were taken at 21 epochs between 6 Aug 1996 and 17 Oct 2011 with, on average, one observation every 259 d. The data were released by Butler et al. (2017), but we used the corrections applied by Tal-Or et al. (2019). We refer to this set of data in the following as season S0, since it does not overlap with our HARPS-N and CARMENES observations and since it is not visibly separated into the seasons of observability. The data show small RV scatter and uncertainty. We considered for the RVs of all three instruments instrumental RV drifts, a barycentric correction following Wright & Eastman (2014), and the small secular acceleration of about 0.12 m s⁻¹yr⁻¹ of Gl 49 (see Table 1).

Fig. 1 RV time series (top panels), GLS periodograms (middle panels), and window functions (bottom panels) of the HIRES (left column), HARPS-N (middle column) and CARMENES (right column) RV data. The different observational seasons are labeled from S0 to S6, and the unusual season S1 is colored in blue. In the periodograms, we mark the periods at 9.37 and 18.86 d with orange vertical lines, the period at 13.85 d with a red vertical line, and the 0.1, 1, and 10% analytical FAP levels by red horizontal dashed lines. The blue vertical dashed line for the HIRES data shows the Nyquist frequency at (159 d)−1 = 0.00629 d−1.

3 Stellar activity

The stellar contribution to RV variations can be very complex and mimic or mask a planetary companion. In order to understand those effectsmore profoundly and to be able to include this knowledge in our RV data modeling, we analyzed the activity level of Gl 49 and collected and derived various additional time series, which are generally understood as tracers for stellar activity effects of the star and on its surface.

3.1 Activity level

Gl 49 showsa moderate activity level, which is expressed by its calcium and Hα indices as shown in Table 1 with = −4.8 dex and pEW(Hα) = −0.04 Å. In the recent ROSAT All-Sky Survey (RASS) source catalog by Boller et al. (2016), the source 2RXS J010318.3+622140 is close to the position of the companion Gl 51. Its count rate of 0.24 ct s⁻¹ leads to an X-ray luminosity of about 1.7 × 10²⁸ erg s⁻¹ for the whole system. This value refers predominantly to Gl 51, however, a close inspection of the images shows an extended source with an elongation toward Gl 49. The observation of Gl 51 (Obs-Id 0742230501) in the XMM-Newton archive contains Gl 49 in its field-of-view. The two sources are clearly separated and the X-ray flux and luminosity can be determined separately. For Gl 49, González-Álvarez et al. (2019) estimate an X-ray luminosity of about 3.9 × 10²⁷ erg s⁻¹ (see Table 1), which accounts for about 20% of the total flux, essentially consistent with the RASS measurements. We then find Gl 49 located in the upper quartile of the early M-dwarf X-ray luminosity distribution function (Schmitt et al. 1995).

3.2 Photometric data

We obtained photometric observations in the framework of the HADES and CARMENES consortia. The basic properties and the temporal distributions of the data are shown in Table 3 and in Fig. A.1.

Simultaneously to the spectroscopic observations during seasons S5 and S6, we collected data within the EXORAP (EXOplanetary systems Robotic APT2 Photometry) program. It is carried out at INAF-Catania Astrophysical Observatory with an 80-cm f/8 Ritchey-Chrétien robotic telescope (Automated Photoelectric Telescope, APT2) located at Serra la Nave on Mt. Etna. BV RI photometry of the star was collected over 157 nights between 7 Jun 2016 and 18 Oct 2017. It covers 499 d with an average of one observation every 2.0–2.3 d for BV RI filters. To obtain differential photometry, we started with an ensemble of about ten stars, the brightest that were close to Gl 49. We checked the variability of each of them by building their differential light curves using the rest of the sample as a reference. In that way, we selected the four least variable stars of the sample for filters B, V, and I (five stars in R). The average rms of the ensemble stars was 15, 15, 24, 21 mmag in the B, V, R, and I filters, respectively. We obtained 250, 239, 230, and 223 data points, respectively, and rejected very few of them by a 5σ clip. We calculated average photometric uncertainties as sky plus Poisson noise. As visible in Fig. A.1, Gl 49 is slightly fainter in S5 than in S6 for EXORAP filter B, but the relationship seems to shift going to longer wavelengths until the star is brighter in S5 for filter I.

We also used the 90 cm Ritchey-Chrétien T90 telescope at Sierra Nevada Observatory (SNO), Spain. It is equipped with a 2k × 2k CCD camera and a field of view of 13.2 × 13.2 arcmin² (Rodríguez et al. 2010). The observations were collected in Johnson V and R filters on 44 nights in 2018 during the period from 11 Jun to 21 Sep. The observations were carried out after the RV campaigns in order to search for a possible planetary transit. The measurements were obtained by the method of synthetic aperture photometry using no binning and in average 25 observations of 10–25 s per night. Each CCD frame was corrected in a standard way for bias and flat-fielding. Different aperture sizes were also tested in order to choose the best one for our observations (16 pixels). A number of carefully selected nearby and relatively bright stars within the frames were used as reference stars. Outlier points due to bad weather conditions or high airmass were removed and a nightly average applied, since no transit event could be detected and no short-periodic signal was expected. The data of the two filters show a correlation with similarvariations (see Fig A.1), but slightly larger errors and a larger apparent amplitude in R band. As visible in Table 3, the average uncertainties exceed the rms of the data of both filters.

Besides the photometry provided by the HADES and CARMENES consortia, we analyzed 11 years of observations from Las Campanas Observatory, Chile, with the ASAS-3N system (All Sky Automated Survey; Pojmański 1997). The 440 data points cover 434 nights from 3 Jun 2006 to 22 Jun 2017. This is a timespan of 4037 d with on average one observation every 9.2 d. The data also cover season S1, S2, S3, and S5 of our RV observations, but show average uncertainties similar to the overall data scatter. The survey delivers a standard magnitude of V_ASAS = 9.631 ± 0.025 mag using a five-pixel aperture.

We also used MEarth data (Berta et al. 2012) obtained at the Fred Lawrence Whipple Observatory on Mount Hopkins, Arizona, USA, and provided to us in three sets (see also Díez Alonso et al. 2019). Set A includes observationson 7 nights from 13 Oct 2008 to 3 Mar 2010. Set B includes observations on 53 nights from 19 Nov 2010 to 23 Jun 2011, and is split into the two seasons of visibility. Set C covers 79 nights from 22 Oct 2011 to 12 Nov 2015 and overlaps with the S1–S4 seasons of our spectroscopic campains (see Fig A.1). We did not find significant periodicities of less than 1 d in the data and applied nightly binning after a 3σ clipping of the values and the uncertainties. The uncertainties are similar to the scatter.

3.3 Spectroscopic indices

We derive a number of additional properties from the HARPS-N and CARMENES spectra. We do not attempt to mix values from different instruments and reduction pipelines. An overview and a visualization of all the time-series data after a 3σ-clip on measurement uncertainties are given in Table 3 and in the left panels of Figs. A.2 (HARPS-N) and A.3 (CARMENES). The data are given in Table A.1.

We calculate the cross-correlation function (CCF) between each observed spectrum and a binary mask constructed from the average spectrum of all observations of Gl 49. The mask is a selection of lines of the combined spectra of certain depth and width, excluding regions of strong atmospheric influences. The CCF is then a representation of the average line of the spectrum in the velocity space and can be approximated by a Gaussian function. The variations over time of the velocity span of the half-maximum power (FWHM) and the maximum depth or contrast (CON) of this fit should therefore correlate directly with velocity variations across the stellar surface. Dumusque et al. (2014) associate such variations with the bright faculae on the stellar surface. We also calculate the bisector inverse slope (BIS, Queloz et al. 2001), which measures the asymmetry of the CCF. This value gives an insight into the structure of the photosphere and should correlate with differences in temperature and pressure, both connected strongly to the activity level of the star. The HARPS-N time-series data show in its temporal distribution already that (i) a long-term variation of >1500 d is present in FWHM and CON, (ii) the BIS value is constant, and (iii) S1 does not stand out in those indices as in the RVs. The CARMENES time-series, on the other hand, are less variable over the different seasons, but match the long-term trend described by the HARPS-N data.

Many magnetically-sensitive line features are present in the optical spectra of both instruments. They are formed in the hot plasma of the chromosphere and hence vary with the strength of the stellar magnetic field. Those features are connected to elements such as calcium, hydrogen, sodium, helium, and iron and are some of the most important tracers for stellar activity in the literature. The Ca II ion shows the Fraunhofer H & K emission lines at 3934 and 3969 Å in a domain of low S/N in the HARPS-N spectra of our target. We quantify those lines, which measure the lower chromosphere (Gomes da Silva et al. 2011), with the S index (CaHK; Duncan et al. 1991). The infrared triplet of the same ion is located in the CARMENES spectral range at 8498, 8542, and 8662 Å (CaIRT). We measure and sum up their fluxes with respect to the continuum. We additionally measure Hα at 6563 Å, which is absorbed rather in the upper part of the chromosphere (Gomes da Silva et al. 2011). We include the measurements of the Na I D1 and D2 absorption doublet at 5890 and 5896 Å. Díaz et al. (2007) proposed that these lines could be used to follow the chromospheric activity level of very active late-type stars and that they provide good complementary information about the conditions in the middle-to-lower chromosphere (Mauas 2000). The HARPS-N indices show the long-period variation >1500 d already described for the RV and CCF index time-series data. In contrast to the CCF indices, the data points show the different behavior of S1 (marked in green in Fig. A.2), but with a strong consistency between the different indices. The CARMENES datasets follow the long-term trend of the HARPS-N sets nicely with an increasing and decreasing activity level for S5 and S6, respectively.

With the publication of the data reduction pipeline SERVAL, the chromatic index (CRX) was introduced as a new tracer for activity phenomena (Zechmeister et al. 2018). For every spectrum, it measures the slope of the RVs calculated for each different echelle order when they are plotted against wavelength. A non-zero and time-variable CRX is therefore indicative of RV variations due to stellar surface phenomena. As shown in Figs. A.2and A.3 the index follows the long-term trend described by, for example the FWHM or CON indices.

3.4 Periodogram analysis

In order to get a first overview of the periodicities in our data sets, we use the Generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009). As thresholds for a signal to be tentative or significant, we use at this stage the 1 and 0.1% False Alarm Probability (FAP), respectively, calculated analytically with the formula by Horne & Baliunas (1986). If the signal shows FAPs below those thresholds, we subtract it by a sinusoidal best-fit and analyze the residuals in the same manner; this procedure is customarily dubbed prewhitening. We show the periodograms of the activity tracers in the right panels of Figs. A.1 (photometry), A.2 (HARPS-N), and A.3 (CARMENES). The red horizontal lines indicate the 0.1, 1, and 10% FAP level, whereas the two orange vertical lines stand for the periodicities at 9.37 and 18.86 d, and the red vertical line marks the period at 13.85 d. All significant signals of the different datasets are given in Table 3 under P_rot, if the signal is close to the 19-d period, and under P_add, if it is not. We mark the first with *, if the signal is distinguishable from the noise but above the 1% FAP level.

The various photometric datasets reveal the period at 18.9 d very well. The first harmonic, on the other hand, is mostly not seen due to the low amplitude of the signal. For the EXORAP filters, the data show an additional strong periodicity at >220 d. We connect it in part to the time-sampling, but also identify the variation in the data point distribution. The SNO photometry shows also a strong periodicity of >150 d. In the noisy ASAS-3 data, we are able to see two statistically insignificant periodicities at 18.94 d and 177 d. Also from MEarth data we are able to clearly detect the rotational periodicity.

It seems that the HARPS-N S1 data are dominated by a strongly varying activity level of Gl 49. This is seen clearly in the large RV scatter of that season (Fig. 1) and in the varying values for the chromospheric indices (Fig. A.2). Since in all cases this data subset (20% of all data points) only adds noise to the periodograms, we exclude it from this analysis and the periodograms in Fig. A.2. Even then, we identify in only four indices low-significance signals at periodicities of 19 or 9.4 d. All sets show at 400 to 425 d a signal which is probably a combination of the long-term periodicity of >1500 d and the strong yearly alias created by the seasons of visibility as seen in the window function, that is, the Fourier transformation of the time-sampling (see Fig. 1). Given the distribution of all the HARPS-N datasets, we assume the star to be magnetically rather active at S2 and by the end of S5. The minimum activity level at approximately Mar 2015 (BJD = 2 457 100 d) coincides with the minimum scatter in the RV measurements. But such signals of longer periods are not found by the periodogram approach.

The data distribution of the CARMENES indices is, on average, very similar between the S4 and S6 seasons, but indicating a slightly less active state of Gl 49 by the end of S6 than in the previous seasons and showing an upward/downward trend for S5/S6. More clearly than for the HARPS-N indices, all their periodograms show features at around 19, and/or 9.4 d, which might in part be due to the higher observing cadence. The alias at 365.25 d shows at longer periods since it is probably as well mixed with the unresolved long-term trend already described.

3.5 Secondary periodicity

The presence and stellar nature of a 19-d periodicity and its first harmonic in Gl 49 is shown by the periodogram analysis of the different activity indices. It is also reflected in the strong symmetry that the periodograms of especially the CARMENES CaIRT and Hα indices showaround this period (see Fig. A.3 and, for the latter, bottom panel of Fig. 2). Besides the important periodicities at P = 18.86 d (orange) and 13.85 d (red), we show the yearly aliases of the first at frequencies P⁻¹ ± 365.25 d⁻¹ (black dashed lines) and the analytical 0.1, 1, and 10% FAP levels (from top to bottom) with the blue horizontal dashed lines. The symmetry results most prominently in two periodogram peaks above the noise level at 14.1 and 28.3 d (green lines), corresponding to frequencies P⁻¹ ± 56.44 d⁻¹. This is indicative of a 56 d-variation of the amplitude of the 19 d-signal. We show the whole CARMENES dataset in the top panel of Fig. 2 in blue. Since this secondary periodicity is not seen in the periodogram but the value is close to 3P, we fold the data to 56.44 d and show it in the middle panel. There, we observe one very strong variation between phase 0.3 and 0.8, proving this secondary 56.4 d-periodicity of the stellar contribution at S5 and S6. We also include the HARPS-N values, after the fit of an adequate offset, which then fall nicely in our picture.

3.6 Time-series data correlations

We identify the rotational periodicity in all photometric and most of the activity index datasets, as well as in the RVs (see Sect. 4). To be able to study their possible correlations, we calculate the phase shifts with respect to the RV time-series data (Suárez Mascareño et al. 2017b). To achieve this, we apply to all sets the best-fit period of the interval of 18 d < P < 20 d of the RV data of HARPS-N (19.11 d), and CARMENES (18.92 d), respectively (see Table 2). The resulting shift for each set is given in Table 3. We expect only the FWHM and CON valuesto be in phase with the RVs. For all other indices, including the photometry, a phase shift of 45° is theoretically expected for a simple one-starspot model. A closer look at the results reveals a great consistency only for the chromospheric indices with phase shifts between 75 and 80° and for the CRX index with 230–240°. We additionally calculate the phase shifts of the photometric data using the HARPS-N periodicity. The results are pointing to a value of 40–80°. The differences for the SNO photometry could be related to an unknown activity event after S6 season.

We then calculate Spearman and Pearson factors and find highly significant linear correlations only for the chromospheric indices. For CaHK/CaIRT and Hα indices we obtain Spearman factors of 0.88 and 0.95 for HARPS-N and CARMENES data, respectively. For the calcium index with the Na I index it is 0.73 and 0.68, and for Hα with the sodium index, we calculate 0.70 and 0.72. We do not find correlations for any other combination of time-series data which is, in most cases, due to the large phase shifts between them. Linear correlations are also not found for RV/CON of CARMENES or CRX/CON of HARPS-N observations, where small phase differences are measured.

To illustrate this procedure, we show the correlations of the RV data with the chromospheric activity indices in Fig. 3. We show, from left to right, the calcium, Hα, and NaI indices for HARPS-N (top panels), and CARMENES (bottom panels) instruments, respectively.

Fig. 2 Top panel: Hα time-series data of CARMENES (blue) and HARPS-N (red) spectra of seasons S4, S5, and S6. Middle panel: data points phase folded to 56.44 d. Bottom panel: periodogram of the CARMENES index as in Fig. A.3. We mark the periods at 18.86 (orange vertical line) and 13.85 d (red vertical line), their yearly aliases at frequencies ν = P−1 ± 365.25−1 d−1 (black vertical dashed lines), and the frequencies at 18.86−1 ± 56.44−1 d−1 (green vertical lines), where the symmetry of the index data strongly imprints. The 0.1, 1, and 10% FAP levels are shown by the blue horizontal dashed lines (from top to bottom).

3.7 Modeling correlated noise

The phase shifts of the activity indices imply that the information about the stellar contribution monitored by each index cannot be implemented easily in the modeling of our RV data. We therefore model the data as correlated noise using GP regression (Rajpaul et al. 2015; Roberts et al. 2013). The covariance matrix of the GP contains all the information we want to apply regarding the correlation of the data points, which is consistent with the picture of active regions rotating with the star and evolving in time. In our case this is a quasi-periodic kernel, similar to a damped oscillator, consisting of a periodic component (described by hyperparameters period P_QP and amplitude K_QP) connected to the rotational period of the star, and a damping factor λ (connected to the evolutionary timescale of surface phenomena). An additional scaling factor w (generally <1) is used to measure the relation of both kernel parts. A small w shows large influence of P_QP and small influence of the exponential decay. Details on this method using the george code by D. Foreman-Mackey are described in Ambikasaran et al. (2015). In another approach, we use the celerite code by the same authors calculating the covariance matrix by a stochastically-driven simple harmonic oscillator (Foreman-Mackey et al. 2017). In this code, we model the rotational periodicity with the two hyperparameters P_SHO and K_SHO, and a hyperparameter P_life, which is supposed to be connected to the lifetime of the active regions.

To compare different models consistently, we use the log-likelihood (ln L) statistics (e.g., Anglada-Escudé et al. 2013). The likelihood function shows the probability of the data matching a certain model, and generalizes the χ² statistics introducing an additional jitter term σ. Details on the definition for a given dataset and noise models are given in, for example Ribas et al. (2018). At this point, we consider a model compared to a second one (including a null model) significant if the difference in ln L exceeds a value of 15, which roughly corresponds to a FAP of 0.1% in our measurements. A difference of Δ ln L < 5 is considered noise. For a detailed discussion on the usage of this statistics for model comparison see Baluev (2013), and for an example implementation of this procedure see Ribas et al. (2018). We provide later (see Sect. 4.2) an empirical determination of the Δ ln L threshold required to claim a confident detection for a certain model.

We further use the emcee code by D. Foreman-Mackey to explore the parameter space with a Markov chain Monte Carlo procedure (MCMC, Foreman-Mackey et al. 2013). For the index time-series data, we apply 10 000 steps on each of 100 walkers. As a comparison, we maximize the ln L using only offset and additional jitter (null-model). In a second step, we then apply the two GP algorithms and their respective covariance matrices. Since we already know the rotational period of Gl 49 and have analyzed both the periodograms and the temporal distribution of the different time-series data, we set boundaries to the uniform priors of the hyperparameters P_QP and P_SHO from 5 to 30 d. To all other priors we apply reasonable boundaries (see Table 5). We show the results of this procedure in Table 3 only if Δ ln L >10 compared to the null model, to include also tentative solutions.

The CaHK/CaIRT and Hα indices and the EXORAP photometry show the most significant and consistent results and retrieve the rotational period very well. But the interpretation of the hyperparameters P_life, λ, and w is more complicated. As already mentioned, w ≫ 1 (e.g., HARPS-N BIS) shows the dominance of the exponential decay over the rotation. On the other hand, λ and P_life are generally not constrained consistently. The latter shows best solutions in the interval 2P_rot–4P_rot for nearly all the indices, confirming the value of 56.4 d, which we found as secondary period in some activity indices of seasons S5 and S6. λ could instead be summarized as an overall larger value working similar for the EXORAP photometry (2P_rot–4P_rot), but showing instead 150–250 d for the other activity indices. The large values for the SNO photometry are due to the stability of the 19 d-variation over the short time-period of the observations.

Fig. 3 Correlations of the HARPS-N (top) and CARMENES (bottom) chromospheric line indices with the RV data. The red ellipses show the theoretical phase-shifted correlation using the phase shifts of periodicities of 19.11 d (HARPS-N) and 18.92 d (CARMENES), respectively, and the rms of the data as the amplitude of the signal. The blue lines represent the best linear fit, i.e., the Spearman factor. At the top of each image, we show the Spearman (SP) and Pearson (PE) factors together with their respective p-values in parenthesis.

4 Radial velocities

4.1 First analysis of individual datasets

To get an overview of the periodicities present in the different RV datasets, we analyze them individually with the GLS method and show the results in Fig. 1 for HIRES, HARPS-N, and CARMENES. The 21 HIRES S0 RV measurements are not able to sample the periodicities that we are interested in, since a lower limit is given by the Nyquist frequency with approximately 160 d. If we search for periodicities down to 1 d, we find, however, a tentative signal at 17.3 d (see Col. 7 of Table 2), which we assume to be connected to the rotational period. In the prewhitening process using the 137 epochs of HARPS-N S1 to S6 observations, we find three strong signals in the range of the rotational period at 18.3, 18.9, and 19.1 d. Furthermore, we find a periodicity close to the first harmonic of the rotation at approximately 9.6 d. The facts that those periodicities are not easily cleaned with simple sinusoidal fits, that they influence each other, and that the supposed harmonic is so prominent point to a rather complex RV contribution from the stellar rotation. The additional significant signal at 13.4 d will be shown to be connected to the orbital period of a proposed planet. If we exclude the S1 dataset in our prewhitening periodogram analysis, the main periodicity of the remaining dataset is now the 9.4 d-signal, as shown in Table 2, and a signal at 14.4 d rather than 13.4 d. The S1 season with the large data scatter shows significant peaks in all three interesting time intervals of the periodogram. The analysis of the 80 CARMENES data points reveals three periodicities at 9.3, 18.9, and 14.4 d. The periodogram is more easily cleaned by sinusoidal fits, indicating a more stable stellar contribution in S4–S6 than in S1–S3. Like in the case of the S2-S6 HARPS-N set, the most prominent periodicity in the CARMENES set is connected to the first harmonic of the rotational period, that is 9.3 d, and a remaining signal at 14.4 d is visible.

4.2 Methodology

Using the knowledge of the periodicities in our datasets, we search in the following for the best model to fit to the RV data, in order to estimate the best parameters for a possible planetary companion. We exclude the first season of the HARPS-N observations, since we assume it to be affected by strong variations in stellar activity. We exclude further the HIRES data, since the time sampling was very poor. Although the rotational period in the RVs is very stable, these datasets add more noise than signal to our data and might distort any fitted parameter. In the following search for the most probable fit parameters, we then use the HARPS-N S2–S6 and CARMENES datasets, that is 80% of the available RV measurements.

We fit different models to be able to compare their probabilities through their ln L values. We use again the emcee code for the MCMC analysis as outlined in Sect. 3.7 with 100 walkers and 10 000 steps. If not mentioned differently, uniform priors with reasonable boundaries (see Table 5) are used, in order to fully exploit the parameter space and to have a different approach to evaluate the probability and significance of the fitted parameters.

We begin by fitting a constant null model, including offsets and jitters for both datasets, and obtaining as a best-fit ln L = −574.0 (see Table 4). Based on this value, we built a periodogram by considering for every period the Keplerian fit (f(P,K,e sin ω, e cos ω, T_per), with P as period, K as amplitude, e as eccentricity, ω as true anomaly, and T_per as time of periastron of the respective orbit) with the largest ln L. As a significance test, we apply bootstrap randomization (Murdoch et al. 1993) using 10 000 permutations of the data points and reach a 0.1 and 1% FAP with Δ ln L = 16.11 and 14.24, respectively. Those numbers confirm the Δ ln L = 15 mentioned above, but are connected to the one-Keplerian model. Nevertheless, in the following, a more detailed model, that is including a correlated noise term (GP) or an additional Keplerian (Kep), is considered tentative or significant if it shows Δ ln L > 14.24 and 16.11, respectively, compared to the former model. We also confirm the noise level of Δ ln L = 5.

Fig. 4 Region of interest of the Δ ln L periodograms of the combined RV datasets of Gl 49 (HARPS-N S2-S6 and CARMENES). Top panel: difference in ln L from a linear null model including additional jitters and offsets, to a model including one additional Keplerian (top panel). The following panel then shows the periodogram of the data with this Keplerian subtracted. From the second panel, from top to bottom: we have subtracted a best-fit Keplerian curve with periods of 9.37, 19.09, 13.85 and 18.91 d, as shown by the broad red vertical lines. The dashed horizontal red lines indicate the 0.1 and 1% FAP values calculated by bootstrap randomization using a single Keplerian model on the datasets (see Sect. 4.2)

4.3 Beat frequency

Disregarding the knowledge of the nature of the RV signals, we fit up to five Keplerians to the datasets. Each of those curves implements five fit parameters, which add to the four parameters already in place for RV offsets and additional jitters (see N_param in Table 4). The increase in ln L is significant until the fourth Keplerian (ln L_4Kep − ln L_3Kep >16.11), where 24 parameters are fitted. The periodograms of this procedure are shown in Fig. 4 for the whole dataset, and the residuals after subtracting the best-fit of one, two and three Keplerians, respectively.

As a result, the four-Keplerian model (see top panel of Fig. 5) is able to clean the periodogram significantly, and it describes the possible stellar contribution with three Keplerian curves (middle panel) with periods of 9.37, 19.09, and 18.91 d. For the sakes of this exercise, the best fit and the largest ln L, this procedure should be as valid as the usual fit of simple- or double-sinusoidal models for the complex and noisy stellar contribution. The remaining periodicity at 13.85 d is hidden behind those three strong signals but identified significantlyas the third Keplerian. In the best MCMC solution, this periodicity shows a RV amplitude of 2.71 m s⁻¹, and an eccentricity of 0.36. Actually, the results are very similar to the outcome of the prewhitening process of the individual datasets using sinusoids in Sect. 4.1. We see strong yearly aliases on all peaks in the periodograms showing the response of the time-sampling as seen in the window function.

Whereas we can explain the signal at 9.4 d as the first harmonic of the rotational period, we show the combinationof the two periods around 19 d in the bottom panel of Fig. 5. The two periodicities are effectively responsible for an amplitude change of the 19 d-signal, which occurs with the beat frequency at ((19.09)⁻¹+(18.91)⁻¹⁾⁻¹ = 2042 d. We identifythis periodicity to be the same as the long-term variation of >1500 d already described, especially in the HARPS-N datasets (see Figs. 1 and A.2).

Fig. 5 Four Keplerian best-fit models including all signals (periods of 9.4, 13.9, 18.9, and 19.1 d; top panel), the three signals connected to the stellar rotation (9.4, 18.9, and 19.1 d; middle panel), and the two signals around 19 d (bottom panel), resulting in the amplitude variation following the beat frequency of 2042 d in the observed time span. The red dots are the HARPS-N S2 to S6 and CARMENES data points.

4.4 Models including correlated RV noise

A more physical approach to model the stellar contribution is the application of a GP noise term. Here, as already explained, we model the stellar contribution as correlated noise through a covariance matrix representing a damped or stochastically driven harmonic oscillator. In Table 4 we can see that the models are significant including also one Keplerian (ln L_1Kep+GP − ln L_GP >16.11) and fitting thereby 13 and 12 parameters for george and celerite codes, respectively. Including a second Keplerian is statistically not important with ln L_2Kep+GP − ln L_1Kep+GP <14.24. With the wide prior boundaries as shown in Table 5, the best model captures both the rotational and the 13.85 d-signal as expected, indicating the difference in stability of those two periodicities. The more flexible quasi-periodic kernel delivers a larger best ln L value in the MCMC analysis than the simple harmonic oscillator kernel, and we reach up to −491.8, resulting in a high significance of Δ ln L = 59.1 for our proposed system and planet. If we compare the results with the model of the four Keplerians of Sect. 4.3, the amplitude of the 13.85 d-signal decreases slightly to 2.54 m s⁻¹. Such a suppression of the RV amplitude is suspected to be a disadvantage of the GP regression (Feng et al. 2016; Ribas et al. 2018).

We carry out some additional tests with the selected best model of one Keplerian plus a GP noise term using the quasi-periodic kernel (see Table 4):

• We assume seasonal RV shifts for each instrument and include in our models RV offsets for every season and instrument. The period ranges that we are interested in should not be affected. As we increase the number of parameters from 13 to 20, we achieve practically no change in ln L. This is certainly partially due to the absorption of the effect by the GP.

• If we apply our model with a circular rather than an eccentric Keplerian orbit, we reach ln L = −496.9. Comparing this circular model to the eccentric model, we reach a difference of Δ ln L = 5.1, which is close to the noise level of Δ ln L = 5.

• If we include HARPS-N S1 data, the RV amplitude drops and the eccentricity slightly rises. We reach ln L = −570.3 for 217 data points, which translates to ln L = −499.3 for 190 data points, indicating a worsening of the fit.

• If we include both HARPS-N S1 and HIRES data, the trends for RV amplitude and eccentricity continue. We reach a translated ln L value of−502.3.

The results of those exercises reflect the stability of both HARPS-N and CARMENES instruments, the unusual behavior of HARPS-N S1 data, and the poor sampling of the HIRES data.

4.5 Final model

To find the most probable values and uncertainties for the 13 free parameters of our selected best model including one Keplerian and a GP noise term, we explore 250 × 50 000 solutions of our MCMC analysis. We use ln L = −519.9 of the GP_QP-term best-fit (see Table 4) as a reference value. We have 10 076 912 solutions above that limit. Including a Keplerian, we consider every solution as tentative which adds to this number at least the 1% FAP of Δ ln L, which is 14.24. The application of such a limit of ln L > −505.7 leaves us with 8 156 650 solutions, which we use for our statistical approach. Our most probable model parameters and their uncertainties are then the median and the rms of the selected solutions, respectively. We show the distribution of all tentative solutions for the fitted parameters in the corner plot in Fig. A.4 and give the final values in Table 5. Additionally, we show the ln L distribution per parameter of interest in Fig. 7. The time series of the final model is shown in Fig. 8 and the phase-folded curve excluding the stellar GP contribution in Fig. 6.

The periods of both the GP term and the Keplerian signal are constrained very precisely, although we can see still some structure well below any noise level in the ln L distribution of the former (see Fig. 7). This is also visible in the corner plot (Fig. A.4), where this parameter shows a slight correlation with λ, which also shows this insignificant bi-modal behavior. The eccentricity of the Keplerian orbit is fitted by e sin ω and e cos ω and adds up to 0.36. We see evidence at this point for the parameters to be different from 0 within Δ ln L < 5 from the peak value. The variable GP-term contribution, which corresponds to the more significant rotational signal, shows a larger impact on the RV data with an amplitude of 3.02 m s⁻¹ than the stable contribution of the Keplerian with 2.52 m s⁻¹. The important values of λ and w are rather poorly constrained. We calculate for the first a value between 100 and 240 d. The median value of w is 4.3. Although beyond the normal limits used for this parameter, it illustrates the importance of the exponential decay, that is the variation in amplitude, of the signal in comparison to its strong periodic behavior. The additional jitters of our best-fit model show a consistent picture and correlate with the uncertainties of the measurements already in place.

Fig. 6 RV data from HARPS-N (green) and CARMENES (red) of Gl 49 of seasons S2 to S6 without the stellar contribution as calculated by the GP term of the final model. The data are folded in phase to 13.85 d to illustrate the impact of the Keplerian orbit as shown by the black line.

Fig. 7 Δ ln L distribution of the different parameters of interest following the MCMC procedure. The zero level is set at the best-fit ln L of the model including only the GPQP term (see Table 4). We show the 0.1 and 1% FAP level calculated by bootstrap randomization using a single Keplerian model (Sect. 4.2) by the black horizontal lines. The black vertical lines in each image indicate the median value and the rms of each parameter considering the solutions below the 1% FAP level. Whereas the green dots indicate all solutions, we plot in red every 1000, and in yellow every 10 000 solutions, in order to get an insight on the density distribution. We note, that for M sin i and a, the uncertainty of the mass of Gl 49 is considered.

5 Discussion

5.1 The planet Gl 49b

We find evidence of a 13.85 d-signal in our RV data. We have shown that the signal and its yearly aliases are not present in any activity index, and that it is not introduced by the time-sampling or the data treatment using a variety of datasets and models. In comparison to the 18.86 d-signal also present in the data, it shows great stability in period, phase, and amplitude over the six years of observations used for the final fit. We therefore assume this periodicity to be caused by the orbital motion of a planetary companion Gl 49b. Its impact on the RV data of its host star is deeply hidden inside the complex stellar activity-induced variations. But with our best noise model, by which we accurately correct for the stellar contribution, the planetary signal is strong, with a significance of Δ ln L = 59.1, and the parameters are very well constrained. We find the companion to have a minimum mass of 5.6 ± 0.7 M_⊕. The super-earth orbits its star in 13.851 ± 0.005 d in an eccentric orbit (e = 0.36 ± 0.10) with a semi-major axis of 0.090 ± 0.001 au. The amplitude it introduces in our RV time-series data is 2.5 ± 0.3 m s⁻¹.

Although not significantly distinguishable from the noise in our GP approach (see Table 4), we consider the eccentricity of the orbit of Gl 49b in our final fit. This is because the distributions of the fitted parameters e sin ω and e cos ω (see Fig. A.4), and of the derived eccentricity (see Fig. 7) are different from zero at the Δ ln L = 5 level. Another argument is that in the prewhitening process of the individual datasets (see Table 2) we detect the yearly aliases at 13.34 and 14.40 d rather than the orbital period of Gl 49b at 13.85 d. Additionally, the eccentricity of the Keplerian signal is very persistent in the variety of models fitted to the whole dataset and to all subsets.

Fig. 8 RV data of Gl 49 of season S2 to S6 of HARPS-N (green) and CARMENES (red), and the final model (black line) and its 1σ uncertainties (gray lines) including a GP term and a quasi-periodic kernel for the stellar contribution and a Keplerian curve. Bottom panel: residuals are shown.

5.2 Stellar rotation

We confirm the rotational period of Gl 49 to be 18.86 d and refine thereby the previously reported values by Suárez Mascareño et al. (2018) and Díez Alonso et al. (2019). The signal induced into the RV data is very stable in phase and period in all the different data subsets, but not in amplitude. The idea is that those variations are caused by surface phenomena, which co-rotate independently of their latitude with the rotational period but appear always at specific longitudes in order to explain the stable phase of the signal (see e.g., Strassmeier 2009). The strong presence of the first harmonic at 9.4 d and its persistence in the prewhitening process shows the non-Keplerian nature of this contribution. The effect can be best explained if we assume the stellar surface to be dark-spot dominated and that the bright hot faculae with their effective convective blueshift have only little influence on the measured RVs. This is expected for low-mass stars and should be also in part responsible for the CCF indices not showing great consistency in this study. In the S1 season, on the other hand, where the activity contribution varies, and where the rotational influence dominates, the impact of faculae is supposed to be stronger.

The periodicity is also detected strongly in the photometric data, where the difference in impact on the different filters shows its non-planetary nature. In general, the photometric data are shifted in phase with respect to the RV data by 40–80°. We detect the rotational signal also in the prewhitening of the CARMENES activity indices, and in four out of seven HARPS-N indices. A closer look with our GP noise model then reveals the presence of the periodicity in all the chromospheric indices. Those datasets show consistent phase shifts of 75–80°, which is in an expected range following Suárez Mascareño et al. (2017b), and slightly smaller than the 120° that we measured for the less active GJ 3942 (Perger et al. 2017b). The CRX index shows a phase shift of 230–240°.

5.3 Evolutionary timescale of dark spots

Signals of a few rotational periods to some hundred d are detected in many of the datasets, which should be at least in part due to the evolutionary time scales of activity phenomena on the surface of Gl 49. Most consistent values of the prewhitening process are delivered by the chromospheric activity indices with periodicities of 400–550 d. Such signals are not found in the RVs. To find better constraints on those timescales we fitted our GP model including P_life and λ, but we still lack a definitive translation of those hyper-parameters to the physical world. For the first, we find values from 40 to 80 d in photometry and chromospheric indices, which match the value of 56.4 d that we found for the amplitude variation of the 19 d-signal in S5 and S6. Whereas the EXORAP photometry finds similar values for λ, the line featuresset the parameter in a range from 150 to 250 d. Since the covariance matrix represents a periodic model with changing apparent amplitude, which is exactly what we see in the RV data, we can most certainly trust the λ = 150 d found by the analysis of those datasets, but scale it down by the Euler number e = 2.71828…, in order to match P_life. Those values then are also in agreement with the decay times found for M dwarfs in Giles et al. (2017). The stellar origin of those periodicitiescan also be seen by the impact shown on the SNO and EXORAP data on this time scale. The contrast of those spots against the photosphere seems to decrease with increasing wavelengths.

5.4 Evolution of stellar activity

The impact of stellar activity on our datasets goes even further, since the RV data shows a varying scatter over the different seasons, which we find to be connected to the amplitude of the stellar signal changing on time scales of at least 1500 d.We can trace this pattern in all activity indices showing that actual evolution of stellar activity is in progress, and the long-term variation is not found in the frequentist analysis of our data. Instead, this signal can be mathematically described by two Keplerian curves of slightly different periods, resulting in a signal with an amplitude that changes on a 2000 d-time scale as a consequence of the beat frequency. As already mentioned, this is exactly what the covariance matrices of the GP noise terms represent and how we explain the high value of w in our GP approach with the quasi-periodic kernel. As a consequence of this study, we assume Gl 49 in S1 to be at a moderate, yet strongly varying activity level, perhaps influenced by bright faculae. The activity level increased until the peak level in S2, which was manifested by an increase in spot number or in the motion of the spots to higher latitudes. Then, until March 2015 (BJD = 2 457 100 d, S3-S4), the activity level dropped. A small increase until S5-S6 was followed by a decreasing level in 2018. This also explains that the activity and planetary signals in all CARMENES data are more stable than in the observations before S4.

6 Conclusions and final remarks

In this study, we analyzed 21.5 yr of time-series data of the M1.5 V-type star Gl 49. Those include 137 epochs of spectroscopic observations with the HARPS-N instrument within the HADES framework, 80 epochs within the M-dwarf program of the CARMENES spectrograph, and additional 21 epochs from archived spectra observed with the HIRES instrument. Only the combination of HARPS-N and CARMENES data and their in-depth analysis made it possible to disentangle the contribution of the exoplanet Gl 49b from the complex activity contributions of its host star. We furthermore studied simultaneous photometric multiband observations made within the EXORAP program, with the T90 telescope at SNO, and using archival MEarth and ASAS photometry.

We discover the super-Earth Gl 49b, which orbits its host star with a period of 13.9 d in an eccentric orbit (e = 0.36) and with a semi-major axis of 0.090 au. In all the different datasets that we have analyzed in this study, this periodicity is only present in the RV time series. The rotational velocity that we are able to calculate from the radius and rotational period of the star is in agreement with the upper limit of 2 km s⁻¹ given for v sini in Table 1, pointing to an orbital inclination close to 90° for an orbit close to the equatorial plane of Gl 49. We therefore calculate the geometric transit probability with 2.0% as a lower limit. Such an event would take approximately 2.5 h and show a depth of 1–2.5 mmagfor a rocky/gaseous nature of Gl 49b with R_p = 1.8 and 3 R_⊕. This is challenging to detect with ground-based telescopes, but could be a good target for TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2015), or the upcoming CHEOPS mission (CHaracterising ExOPlanet Satellite; Broeg et al. 2014). We furthermore calculate an equilibrium temperature of approximately 350 K for the planet assuming zero albedo and using the radiative equilibrium temperature T_{p, eq.} = . Gl 49b orbits well interior to the “recent Venus” habitable zone at approximately 0.18 au (see Table 1). In the search for exoplanets around M dwarfs, Gl 49b has characteristics similar to other recently discovered planets with its minimum mass of 5.6 M_⊕ and orbital period of 13.9 d (see Affer et al. 2019).

The radial-velocity signal of the planet is hidden among strong signals of various time scales induced by the rotating host star and its evolving surface elements. Recently, such cases have been discovered quite often in the search for exoplanets orbiting early M-type stars (e.g., Affer et al. 2016; Perger et al. 2017b), which shows the importance of the understanding of the stellar contributionto all different kinds of observed time-series data. Whereas there are no signs of differential rotation on the stellar surface of Gl 49, the activity level varies significantly over the observed timespan. This is reflected in a long-term variation, which we suspect to be connected to the motion and reduction of the number of dark spots. We do not have sufficient data at hand, but if this behavior were periodic as in the magnetic cycle on the Sun, we would estimate its period to be >1500 d, which is consistent with recent findings for early M dwarfs (Suárez Mascareño et al. 2018). We also detect a periodicity of some 40–80 d, which we explain with the timescale for the evolution of spots and faculae. Interestingly, those signals are visible in all our different datasets no matter if they are rather measurements of phenomena produced by actual motions (CCF), or by stellar activity (lines, photometry), which are rather connected to changes in temperature or pressure. From 2012 to 2018, the star went from an unusual season at the beginning up to a peak of activity in 2013–2014. After passing a local minimum in 2015, the level of activity rose to a less-active local maximum in 2016–2017 and is since then again slowly losing strength. In the framework of the CARMENES M-dwarf program, Gl 49 will be observed in the up-coming semesters in order to monitor this interesting behavior in all the different time-series data.

Acknowledgements

M.P., I.R., J.C.M, D.B., 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. G.S. acknowledges financial support from “Accordo ASI-INAF” No. 2013-016-R.0 July 9, 2013 and July 9, 2015. The Italian Telescopio Nazionale Galileo (TNG) is operated on the island of La Palma by the INAF – Fundación Galileo Galilei at the Roque de los Muchachos Observatory of the Instituto de Astrofísica de Canarias (IAC). 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. CARMENES is an instrument for the Centro Astronómico Hispano-Alemán de Calar Alto (CAHA, Almería, Spain). CARMENES is funded by the German Max-Planck-Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF FICTS-2011-02 funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Insitut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the Spanish Ministry of Science through projects RYC2013-14875, AYA2015-69350-C3-2-P, AYA2016-79425-C3-1/2/3-P, ESP2016-80435-C2-1-R, ESP2017-87143-R, ESP2017-87676-C05-1/2/5-R, and AYA2017-86389-P, the German Science Foundation through the Major Research Instrumentation Program and DFG Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. Additional support was provided by the European Union FP7/2007-2013 program under grant agreement No. 313014 (ETAEARTH), the Progetto Premiale INAF 2015 FRONTIER, and the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). This work has made 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.

Appendix A Observational log, data distributions, periodograms, and MCMC solutions

Fig. A.1 Differential photometric data of Gl 49 in mmag (left panels; we note the inverted y axis) and GLS periodograms (right panels) including EXORAP, SNO, ASAS, and MEarth. The vertical orange lines indicate periods of 9.37 and 18.86 d, the red vertical line a period of 13.85 d, and the horizontal dashed red lines the 0.1, 1, and 10% analytical FAP.

Fig. A.2 Seven activity indices calculated from HARPS-N spectra. Left: time series of (from top to bottom) the CaHK, Hα, NaI, FWHM, CON, BIS, and CRX indices. Right: periodograms of those datasets excluding HARPS-N S1 (marked in green on the left). For more details see Fig. A.1.

Fig. A.3 Same as Fig. A.2 but for the CARMENES data.

Fig. A.4 Corner plot of the fitted parameters of the MCMC solutions using a Keplerian orbit (parameters P, K, e sin ω, e cos ω, Tper) and a GP noise term (hyper-parameters PQP, KQP, w, λ) to the RVs of HARPS-N S2 to S6 and CARMENES of Gl 49. The solutions included in this plot exceed the ln L value of a fit with a GP term only (see Table 4) plus the 1% FAP value that we calculate empirically in Sect. 4.2.

References

All Tables

All Figures

Fig. 1 RV time series (top panels), GLS periodograms (middle panels), and window functions (bottom panels) of the HIRES (left column), HARPS-N (middle column) and CARMENES (right column) RV data. The different observational seasons are labeled from S0 to S6, and the unusual season S1 is colored in blue. In the periodograms, we mark the periods at 9.37 and 18.86 d with orange vertical lines, the period at 13.85 d with a red vertical line, and the 0.1, 1, and 10% analytical FAP levels by red horizontal dashed lines. The blue vertical dashed line for the HIRES data shows the Nyquist frequency at (159 d)−1 = 0.00629 d−1.

In the text

Fig. 2 Top panel: Hα time-series data of CARMENES (blue) and HARPS-N (red) spectra of seasons S4, S5, and S6. Middle panel: data points phase folded to 56.44 d. Bottom panel: periodogram of the CARMENES index as in Fig. A.3. We mark the periods at 18.86 (orange vertical line) and 13.85 d (red vertical line), their yearly aliases at frequencies ν = P−1 ± 365.25−1 d−1 (black vertical dashed lines), and the frequencies at 18.86−1 ± 56.44−1 d−1 (green vertical lines), where the symmetry of the index data strongly imprints. The 0.1, 1, and 10% FAP levels are shown by the blue horizontal dashed lines (from top to bottom).

In the text

Fig. 3 Correlations of the HARPS-N (top) and CARMENES (bottom) chromospheric line indices with the RV data. The red ellipses show the theoretical phase-shifted correlation using the phase shifts of periodicities of 19.11 d (HARPS-N) and 18.92 d (CARMENES), respectively, and the rms of the data as the amplitude of the signal. The blue lines represent the best linear fit, i.e., the Spearman factor. At the top of each image, we show the Spearman (SP) and Pearson (PE) factors together with their respective p-values in parenthesis.

In the text

Fig. 4 Region of interest of the Δ ln L periodograms of the combined RV datasets of Gl 49 (HARPS-N S2-S6 and CARMENES). Top panel: difference in ln L from a linear null model including additional jitters and offsets, to a model including one additional Keplerian (top panel). The following panel then shows the periodogram of the data with this Keplerian subtracted. From the second panel, from top to bottom: we have subtracted a best-fit Keplerian curve with periods of 9.37, 19.09, 13.85 and 18.91 d, as shown by the broad red vertical lines. The dashed horizontal red lines indicate the 0.1 and 1% FAP values calculated by bootstrap randomization using a single Keplerian model on the datasets (see Sect. 4.2)

In the text

Fig. 5 Four Keplerian best-fit models including all signals (periods of 9.4, 13.9, 18.9, and 19.1 d; top panel), the three signals connected to the stellar rotation (9.4, 18.9, and 19.1 d; middle panel), and the two signals around 19 d (bottom panel), resulting in the amplitude variation following the beat frequency of 2042 d in the observed time span. The red dots are the HARPS-N S2 to S6 and CARMENES data points.

In the text

	Fig. 6 RV data from HARPS-N (green) and CARMENES (red) of Gl 49 of seasons S2 to S6 without the stellar contribution as calculated by the GP term of the final model. The data are folded in phase to 13.85 d to illustrate the impact of the Keplerian orbit as shown by the black line.
In the text

Fig. 7 Δ ln L distribution of the different parameters of interest following the MCMC procedure. The zero level is set at the best-fit ln L of the model including only the GPQP term (see Table 4). We show the 0.1 and 1% FAP level calculated by bootstrap randomization using a single Keplerian model (Sect. 4.2) by the black horizontal lines. The black vertical lines in each image indicate the median value and the rms of each parameter considering the solutions below the 1% FAP level. Whereas the green dots indicate all solutions, we plot in red every 1000, and in yellow every 10 000 solutions, in order to get an insight on the density distribution. We note, that for M sin i and a, the uncertainty of the mass of Gl 49 is considered.

In the text

	Fig. 8 RV data of Gl 49 of season S2 to S6 of HARPS-N (green) and CARMENES (red), and the final model (black line) and its 1σ uncertainties (gray lines) including a GP term and a quasi-periodic kernel for the stellar contribution and a Keplerian curve. Bottom panel: residuals are shown.
In the text

	Fig. A.1 Differential photometric data of Gl 49 in mmag (left panels; we note the inverted y axis) and GLS periodograms (right panels) including EXORAP, SNO, ASAS, and MEarth. The vertical orange lines indicate periods of 9.37 and 18.86 d, the red vertical line a period of 13.85 d, and the horizontal dashed red lines the 0.1, 1, and 10% analytical FAP.
In the text

	Fig. A.2 Seven activity indices calculated from HARPS-N spectra. Left: time series of (from top to bottom) the CaHK, Hα, NaI, FWHM, CON, BIS, and CRX indices. Right: periodograms of those datasets excluding HARPS-N S1 (marked in green on the left). For more details see Fig. A.1.
In the text

	Fig. A.3 Same as Fig. A.2 but for the CARMENES data.
In the text

Fig. A.4 Corner plot of the fitted parameters of the MCMC solutions using a Keplerian orbit (parameters P, K, e sin ω, e cos ω, Tper) and a GP noise term (hyper-parameters PQP, KQP, w, λ) to the RVs of HARPS-N S2 to S6 and CARMENES of Gl 49. The solutions included in this plot exceed the ln L value of a fit with a GP term only (see Table 4) plus the 1% FAP value that we calculate empirically in Sect. 4.2.

In the text