© ESO 2020

1 Introduction

More than 4200 exoplanets have been confirmed so far¹. A significant fraction have been discovered with the radial velocity (RV) method. In the past decades, the development of new high-precision spectrographs allowed probing a large variety of planets with minimum masses of several Jupiter masses down to only 0.7M_⊕ for YZ Cet b, which is the least massive planet detected so far with the RV technique (Astudillo-Defru et al. 2017; Stock et al. 2020). One such high-precision spectrograph is the CARMENES instrument (Quirrenbach et al. 2014, 2018), which is used to conduct a survey for detecting exoplanets around M dwarfs, which are the most abundant stars of our Galaxy (Kroupa 2001; Chabrier 2003; Henry et al. 2006). The detection of the large number of exoplanets has resulted in the discovery of exotic new types of planets that have no counterpart in our own Solar System, such as super-Earths (M = 1.9–10 M_⊕; Rivera et al. 2005; Valencia et al. 2007; Charbonneau et al. 2009). These super-Earths are abundant around M dwarfs (Dressing & Charbonneau 2015).

The detection of planets close to or inside the habitable zones (HZ, see Kasting et al. 1993; Kopparapu et al. 2013) of their parent stars is of particular interest. With the current technology, M dwarfs are ideal targets for detecting such temperate planets because the HZ of these stars corresponds to a relatively small orbital radius. The lower host star masses result in a higher Doppler amplitude (higher by a few m s⁻¹) than those of more massive stars, which can be measured by current techniques. However, M dwarfs tend to be very active (Johns-Krull & Valenti 1996; Delfosse et al. 1998; Mohanty & Basri 2003; Reiners et al. 2012). The activity can make the detection of small planets difficult by inducing distortions in the shape of the spectral line profiles; thismimicks a planetary signal (Queloz et al. 2001; Desort et al. 2007; Barnes et al. 2011; Robertson et al. 2014, 2015).

Various methods can be used to distinguish stellar astrophysical signals from planet-induced signals. Photometric observations, ideally contemporaneous with the RV observations, as well as different spectral activity indicators can be used to derive more information on the stellar rotation period and activity-induced RV variations. In addition, many novel techniques have been developed to analyze the coherence of a signal, for example, Bayesian-stacked periodograms (Mortier et al. 2015; Mortier & Collier Cameron 2017), growth of the Lomb–Scargle power, or the evolution of the significance (Hatzes 2013; Ribas et al. 2018; Reichert et al. 2019). These tools can provide strong indications that a signal has a nonplanetary origin because an RV signal caused by Keplerian motion should be coherent and long-lived. When both planetary signals and activity contribute significantly to the RV variations, it can be necessary to simultaneously fit for these signals using Gaussian process (GP) regression (Rajpaul et al. 2015) or similar models, such as sinusoids (Boisse et al. 2011; Dumusque et al. 2012). Modeling the stellar activity simultaneously with the Keplerian fit is essential because this contamination can have a significant effect on the derived planetary parameters (see, e.g., Stock et al. 2020).

In the following, we present a detailed analysis of photometric and spectroscopic data of GJ 251, HD 238090, and Lalande 21185. For GJ 251, Butler et al. (2017) reported a possible planet candidate at a period of 1.74 d, but with the more precise CARMENES data, we cannot confirm this claim. Lalande 21185 is the brightest M dwarf in the northern hemisphere, the fourth closest main-sequence star system after α Centauri, Barnard’s star, and CN Leo, and the third closest planetary system. Lalande 21185 has a remarkable history regarding former planet claims. Those by van de Kamp & Lippincott (1951) and Gatewood (1996) were based on astrometric data, but have never been confirmed independently. Later, Butler et al. (2017) reported that Lalande 21185 has a planet candidate with an orbital period of 9.87 d. A recent study by Díaz et al. (2019) was unable to provide evidence for these previous planet claims. However, Díaz et al. (2019) announced the discovery of a super-Earth planet orbiting Lalande 21185 with a period of 12.93 d. The analysis of our CARMENES data agrees with the findings from Díaz et al. (2019) and confirms a single planet orbiting Lalande 21185. HD 238090 has no reported planet to date.

In Sects. 2 and 3, we describe the data, instruments, and methods we used within this study, while in Sect. 4 we compile the basic stellar properties of GJ 251, HD 238090, and Lalande 21185. We then analyze our photometric and RV data for the three stars in Sects. 5, 6, and 7, and provide a star-by-star discussion in Sect. 8 and a general summary in Sect. 9.

2 Data

2.1 High-resolution spectroscopy

CARMENES

GJ 251, HD 238090, and Lalande 21185 were observed as part of our CARMENES² guaranteed-time observation program (GTO) to search for exoplanets around M dwarfs (Reiners et al. 2018a). CARMENES is a double-channel échelle spectrograph installed at the 3.5 m telescope of the Calar Alto Observatory in Almería, Spain. Details regarding the instrument and its performance are given in Quirrenbach et al. (2014, 2018), Reiners et al. (2018a), and Trifonov et al. (2018). The data were processed with the standard pipelines and were reduced with caracal (Caballero et al. 2016b). The RVs obtained with serval (Zechmeister et al. 2018) were corrected for barycentric motion, secular perspective acceleration, instrumental drift, and nightly zero-point variations (Trifonov et al. 2018, 2020; Tal-Or et al. 2019). Table 1 shows a summary of the CARMENES visual arm (VIS) RVs and their overall quality. The median exposure times in the VIS channel were 509, 1000, and 95 s, resulting in a median signal-to-noise (S/N) of 116, 154, and 132 for GJ 251, HD 238090, and Lalande 21185, respectively.

In the CARMENES near-infrared (NIR) data, the scatter was not sufficiently small for the RV analysis of the planetary signals in this work; it was on the order of a few m s⁻¹ (see Bauer et al. 2020, for a detailed analysis of the performance of CARMENES). The RV time series and their uncertainties for the CARMENES VIS data of GJ 251, HD 238090, and Lalande 21185 are listed in Tables D1–D3 together with some activity indicators.

HIRES

The High-Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) is installed at the Keck I telescope in Hawai’i, USA. HIRES uses the iodine cell technique (Butler et al. 1996) to obtain RV measurements with a typical precision of a few m s⁻¹. We used archival HIRES data for GJ 251 and Lalande 21185 to confirm the planetary signals and to extend the time baseline, and to search for long-period signals. For our analysis, we used the HIRES data corrected by Tal-Or et al. (2019), which account for nightly zero-point offsets and an instrumental jump in 2004, which is an improvement over the original data reduction by Butler et al. (2017). Details on the quality of the data are given in Table 1. The median exposure times for GJ 251 and Lalande 21185 were 500 and 135 s, respectively.

SOPHIE

We also used RV data for Lalande 21185 obtained with the SOPHIE instrument (Perruchot et al. 2008). These data were made public by Díaz et al. (2019), and further information on the acquisition and properties of these data is provided in their study. We show a summary of the quality of the RV data in Table 1.

2.2 Photometry

We carried out a contemporaneous photometric follow-up of GJ 251 and HD 238090 during 2018 and 2019. We also compiled photometricdata publicly available as described below.

T90

We monitored GJ 251 and HD 238090 in the Johnson V and R bands with the T90 telescope at the Observatorio de Sierra Nevada (OSN) in Granada, Spain. The T90 telescope is a 90 cm Ritchie–Chrétien telescope equipped with a 2k × 2k pixel VersArray CCD camera, with a field of view of 13.2 × 13.2 arcmin² (Rodríguez et al. 2010). The observations of GJ 251 and HD 238090 were carried out on 42 nights from October 2018 to February 2019 and on 53 nights from February 2019 to July 2019, respectively. The typical number of exposures per night and target was around 35. We did not apply any binning, corrected each CCD frame in a standard way for bias and flat-fielding with IRAF, and selected the best aperture sizes and reference stars for the synthetic aperture photometry. In particular, we used the same aperture size as in Perger et al. (2019).

TJO

Observations of GJ 251 and HD 238090 with the 80 cm Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec in Lleida, Spain, were conducted using a Johnson R filter and its main imaging camera LAIA, a 4k × 4k back-illuminated CCD with a pixel scale of 0.4 arcsec and a field of view of 30 × 30 arcmin². The TJO data for GJ 251 were collected between February and November 2019 during 157 nights and for HD 238090 between February and November 2019 during 149 nights. We obtained several batches of five images per night. The images were calibrated with darks, bias, and flat fields with the icat pipeline (Colome & Ribas 2006). Differential photometry was extracted with AstroImageJ (Collins et al. 2017) using the aperture size and the set of comparison stars that minimized the root mean square (rms) of the photometry. Data with low S/N due to bad weather conditions or high airmass were removed. For GJ 251, we removed 392 low S/N measurements from the initial 2746 data points and for GJ458A, we removed 477 from the initial dataset of 6207 measurements. These correspond to the measurements for which the S/N of the target is below 30% of the best measurement. The resulting light curves were binned to one data point per hour.

LCO

We observed GJ 251 on 44 epochs using the 40 cm telescopes of Las Cumbres Observatory (LCO) in the V band at the Teide, Haleakala, and McDonald observatories between 13 January and 3 March 2019. The telescopes are equipped with a 3k × 2k SBIG CCD camera with a pixel scale of 0.571 arcsec, providing a field of view of 29.2 × 19.5 arcmin². We acquired 50 individual exposures of 30 s per epoch. Weather conditions were mostly clear, and the average seeing varied from 1.5 arcsec to 3.0 arcsec. Raw data were processed using the banzai pipeline (McCully et al. 2018)³, which includes bad pixel, bias, dark, and flat-field corrections for each individual night. We performed aperture photometry for GJ 251 and three reference stars in the field and obtained the relative differential photometry. We adopted an aperture of 13 pixels (7.4 arcsec), which minimized the dispersion of the differential light curve.

TESS

The Transiting Exoplanet Survey Satellite (TESS) is a space-borne instrument that searches for transiting planets around nearby stars (Ricker et al. 2015). The primary mission goal consists of observations of 26 sectors with 24 × 96 deg² in the northern and southern hemisphere, which are still ongoing. Each sector is observed for about 28 d. We obtained for all three targets of this work the pre-search data conditioning simple aperture photometry (PDCSAP) light curves. These are provided by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at the Mikulski Archive for Space Telescopes (MAST)⁴.

MEarth

We used data of HD 238090 from the seventh data release (DR7⁵) of the MEarth project (Berta et al. 2012). The MEarth project is an all-sky transit survey that has been conducted since 2008. It consists of 16 robotic 40 cm telescopes, 8 located in the Northern Hemisphere at the Fred Lawrence Whipple Observatory in Arizona, USA, and the other 8 in the Southern Hemisphere located at Cerro Tololo Inter-American Observatory, Chile. The project monitors several thousand nearby mid- and late-M dwarfs over the whole sky. Each telescope is equipped with a 2k × 2k CCD that provides a field of view of 26 × 26 arcmin². MEarth uses an RG715⁶ long-pass filter, except for the 2010–2011 season, when an I_715−895 interference filter was chosen.

NSVS

The Northern Sky Variability Survey (Woźniak et al. 2004, NSVS) was a robotic survey that primarily targeted the northern sky with telephoto lenses located at the Los Álamos National Laboratory in New Mexico, USA. The survey provided data for 14 million objects in the magnitude range between 8 mag to 15.5 mag. For details on the instrumental setup and the conducted observations, we refer to the survey paper (Woźniak et al. 2004). We used public NSVS data for HD 238090, which we obtained from their public webpage⁷.

SuperWASP

For GJ 251 we used public data processed and collected by the Wide-Angle Search for Planets (WASP) survey (Pollacco et al. 2006)⁸, in particular SuperWASP-North at the Observatorio del Roque de los Muchachos in La Palma, Spain. SuperWASP-North consisted of one wide-field array of eight cameras, each with a 200 mm, f/1.8 lens, a broadband filter spanning the wavelength range between 400 nm and 700 nm, and a 2k × 2k CCD. The resulting plate scale was 13.7 arcsec pixel⁻¹.

3 Methods

3.1 Periodograms

We used generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) to assess significant periodicities in the photometric and spectroscopic data. We applied the normalization as given in Zechmeister & Kürster (2009), which is abbreviated as P_ZK throughout this work. For each periodogram, we computed false-alarm probabilities (FAPs) by applying bootstrapping with n = 10 000 iterations. Our detection threshold for a signal deemed to be significant was at an FAP < 0.001. The uncertainties on the periods of significant GLS signals were estimated from the local χ² curvature by the GLS routine.

To assess the coherence of a periodic signal over the observation time, we used the stacked-Bayesian GLS periodogram (s-BGLS; Mortier et al. 2015; Mortier & Collier Cameron 2017). The Bayesian GLS periodogram allows the comparison of probabilitiesof periodic signals in the data, while the stacking examines the coherence of the signal with an increasing number of observations. As in Mortier & Collier Cameron (2017), we normalized all s-BGLS periodograms to their respective minimum values, which means that the probability of each signal and its growth or decrease over time is a relative measure compared to the lowest probability obtained within one calculated s-BGLS over a specific period range.

3.2 Modeling of RV and photometric data

For the modeling, we used juliet (Espinoza et al. 2019), which allows the fitting of photometric and RV data by searching for the global posterior maximum based on the evaluation of the Bayesian log-evidence () within a provided prior volume of the fitting parameters. juliet allows us to statistically compare models with different numbers of parameters within a Bayesian framework through the log-evidence, which includes the model complexity and the number of degrees of freedom within its assessment. Following Trotta (2008), a model is considered as a significant improvement if . The juliet calculation of the log-evidence is conducted with nested sampling algorithms. In particular, we used the dynamic nested sampling algorithm dynesty (Speagle 2020).

We used radvel (Fulton et al. 2018) to model Keplerian RV signals, and george (Ambikasaran et al. 2015) for GP modeling of both photometric and RV data. In all cases, we used an exp-sin-squared kernel multiplied with a squared-exponential kernel, which is included as a default kernel within juliet. This kernel, also known as the quasi-periodic (QP) kernel, has the form (1)

where σ_GP is the amplitude of the GP component given in parts per million (ppm) for photometric data or m s⁻¹ for RV data, Γ is the amplitude of the GP sine-squared component and is dimensionless, α is the inverse length-scale of the GP exponential component given in d⁻², P_rot the period of the GP QP component given in d, and τ is the time lag. This choice of kernel represents one part of our prior knowledge, as it provides the framework of how an effective model of stellar activity should fit the data. The timescale P_dec in days of the exponential decay can be approximated with (2)

The α_GP parameter is of particular interest with regard to the stability of a QP signal. A smaller α describes a more stable periodic signal in which data points are more strongly correlated with each other. For a review and a detailed description of each kernel hyperparameter and a possible physical interpretation, we refer to Angus et al. (2018).

The evaluation of the GP likelihood with george is computationally expensive and scales as N ln N, where N is the number of data points (Ambikasaran et al. 2015). For the derivation of the stellar rotation, we searched for periods on timescales of days. To do this, it is reasonable to create nightly bins of the photometric data. This reduces the computation time of the GP log-likelihood evaluation and short-term variations based on the jitter of the star.

For the photometric analysis, we applied distinct GP hyperparameters for the amplitudes σ_GP and Γ, to account for the effect that stellar activity depends on wavelength, but we used global GP hyperparameters for the timescale of the amplitude modulation and the rotation period. In addition, we fit an offset and a jitter term (in quadrature to the diagonal of the resulting covariance matrix of the GP) for each data set. Table A.1 shows the priors of the photometric GP analysis.

For the final RV analysis, we applied global GP hyperparameters. A statistical comparison with models using distinct GP hyperparameters for each RV instrument did not show any significant improvement in log-evidence. For each data set, we fit an offset and a jitter term. Our priors for the RV GP analysis are provided in Table A.3.

3.3 De-aliasing

Aliases are spurious signals caused by the sampling of the data, which are often indistinguishable from the true signal. The significance of a signal or the goodness of a fit is not a sufficient criterion to differentiate between true signals and alias signals, especially incases of non-optimal sampling, where the results of these metrics can be similar. For example, it is a common misconception that peaks close to one day are always alias frequencies, but a priori, it is not clear which of the peaks represents the true frequency of the signal and which represents the alias (Dawson & Fabrycky 2010). Alias frequencies can be calculated by f_a = f_t ± mf_s, where f_t is the assumed true frequency, f_s the sampling frequency, and f_a the alias frequency. Because RV measurements are usually taken with a rather irregular sampling (Garcia-Piquer et al. 2017), more than one sampling frequency is often apparent in the window function of the data. This results in several peaks at alias frequencies related to the different sampling frequencies, and can make it even harder to distinguish the true underlying signal.

We used the AliasFinder (Stock & Kemmer 2020)⁹ to confirm that the assumed planetary signal is the true signal and not an alias. The method on which the AliasFinder is based is described in Dawson & Fabrycky (2010), Stock & Kemmer (2020), and Stock et al. (2020). For each frequency under consideration, AliasFinder simulates 1000 data sets based on the true sampling of the observed data and inserts one sinusoidal signal with one of the frequencies. AliasFinder also includes a noise contribution based on the RV jitter of the star, which we made use of for our analyses in this paper. We compared the resulting ensemble periodograms for each simulated frequency to the periodogram obtained from the observed data. Peak position, power, and phase are the parameters compared in this test. If the ensemble of the periodograms of one simulated frequency reproduces the data periodogram significantly better than the periodograms of the other simulated frequencies, then the most probable planetary period has successfully been identified.

4 Stellar properties

Photospheric parameters, such as effective temperature, surface gravity, and metallicity, were determined by Schweitzer et al. (2019) by fitting an updated set of PHOENIX-ACES atmosphere models (Husser et al. 2013) to high-resolution CARMENES spectra. These updated PHOENIX models incorporated the latest solar abundances, molecular and atomic line lists, and a new equation of state (Meyer 2017), which were especially designed to treat low-temperature stellar atmospheres. The parameters were determined assuming a rotational velocity of v sini = 2 km s⁻¹ (Reiners et al. 2018b). To reduce degeneracies between the parameters, Schweitzer et al. (2019) constrained the surface gravity logg with the help of evolutionary models (PARSEC, Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014) and stellar ages estimated by Passegger et al. (2019). The actual ages of the three investigated stars are probably older than tabulated, as derived from a new kinematics analysis (Cortés-Contreras et al., in prep.). The galactocentric space velocities in Table 2 were calculated from the latest Hipparcos and Gaia DR2 proper motions and parallaxes (van Leeuwen 2007; Gaia Collaboration 2018) and absolute RVs of Lafarga et al. (2020), following the approach of Montes et al. (2001) and Cortés-Contreras (2016).

Physical parameters, such as luminosity, radius, and mass, were derived by Schweitzer et al. (2019). Cifuentes et al. (2020) exhaustively described the luminosity determination in M dwarfs. The radius (and hence mass) for Lalande 21185 was an outlier in Schweitzer et al. (2019) because the photometry for this star was of low quality, suggesting an uncertain and too low luminosity. They derived a stellar radius and mass of 0.3587 ± 0.0157 R_⊙ and 0.355 ± 0.019 M_⊙, respectively.We therefore used a slightly different approach to derive the mass and radius. We used its radius 0.3921 ± 0.0037 R_⊙, which was derived by Boyajian et al. (2012) using the interferometric angular diameter. Applying the same empirical mass-radius relationship as was used for the other two targets, we derived a stellar mass of 0.390 ± 0.011 M_⊙, which agrees better with the typical parameters of the ensemble. The detailed stellar parameters of all three stars and their references are given in Table 2.

The M0.0 V star HD 238090 is the primary component of a wide binary. The secondary component is the M3.0 V star GJ 458 B, with a stellar mass of 0.230 ± 0.005M_⊙. The angular separation of the two components of 14.68 ± 0.44 arcsec (Cortés-Contreras 2016) results in a projected separation of approximately 224 au. We computed the stellar mass of the secondary component using the mass-luminosity-metallicity relation of Mann et al. (2019). Based on the masses and projected minimum separation, we estimated the minimum orbital period of this binary to be longer than 3700 yr. This long binary period agrees with the 34 observations between 1955 and 2015 tabulated in the Washington Double Star Catalog (Mason et al. 2001), which do not indicate any change in the position angle.

5 GJ 251

5.1 Photometric monitoring

For a significant fraction of the CARMENES RV observations of GJ 251, we obtained quasi-simultaneous photometry with the T90, TJO, and LCO telescopes. We combined these data with public data from SuperWASP. A joint GLS periodogram analysis, where we fitfor offsets and jitter of each data set, indicated significant signals at periods of 30 d, 70 d, and 120 d. However, a sinusoidal model, as used in the GLS analysis, is an imperfect description of stellar activity, which is often better represented by a QP signal. Therefore we fit a more sophisticated model to the photometric data in the form of a QP GP to derive the stellar rotation period.

Our GP analysis of the photometry of GJ 251 based on our T90, TJO, LCO, and SuperWASP data resulted in a bimodal distribution for the rotational period with posterior solutions around 120 and 60 d. We plot the informative GP α-period diagram (α_GP versus P_GP) in Fig. 1. This plane of parameters shows the decay-timescale over the rotation period, and it is useful for identifying whether stronger correlated noise (small α) favors a certain periodicity (see also Stock et al. 2020 for a more detailed explanation). Within this plane, we identified that the likelihood and number of posterior samples at 120 d is higher than that of the posterior samples around 60 d. Furthermore, the α values of the 120 d signal converge toward our prior boundary of 10⁻¹⁰ d⁻², representinga decay timescale longer than 70 000 d, which indicates a stable periodic signal over the entire time of observations. We determined the rotational period for each posterior solution and derived d and d.

The latter we regard formally as the derived rotational period of GJ 251 because on average, its likelihood of posterior samples is higher than the former solution, because of the stronger coherence of the signal, and because the 60 d signal can be explained as the first harmonic of a signal with a fundamental period of about 120 d. A stronger coherence of signals related to stellar activity would be expected for M dwarfs because the spot lifetime increases with decreasing effective temperature (Giles et al. 2017; Shapiro et al. 2020). Additionally, if GJ 251 is a slowly rotating star, which means that it is relatively inactive, there is evidence that faculae, which are in general longer-lived then starspots, are dominant surface features (see Shapiro et al. 2020, and references there). These long-lived faculae, in particular, affect the photometric variability of the star (Reinhold et al. 2019) and less so the RVs, which are typically spot-dominated. For this reason, among others, the same decay timescales α_GP of the rotational signal between RV and photometric data should not be assumed. We show the binned photometric data overplotted with the median GP model and its uncertainties in Fig. 2.

Based on our estimate of the stellar rotation period, we used Eqs. (1) and (2) by Suárez Mascareño et al. (2018) to estimate and based on this, the expected RV semiamplitude of the stellar rotational signal. By propagating all uncertainties of the parameters given by Suárez Mascareño et al. (2018) for M0-M3 stars and our measurement uncertainty of the rotation period (in the form of their actual distributions), we derived for the medians and 1σ uncertainties (mean at − 5.83) and (mean at 3.59 m s⁻¹).

Fig. 1 Posterior distribution in the αGP vs. Prot plane of the GP fit to the combined photometric data of T90, TJO, LCO, and SuperWASP for GJ 251. The color-coding shows the log-likelihood normalized to the highest achieved log-likelihood value within the posterior sample. Gray samples indicate solutions with Δ ln L lower than 10.

5.2 Spectroscopic activity indicators

We analyzed a number of spectral activity indicators for GJ 251 obtained from the CARMENES spectra using the indicators provided by serval, which includes the chromatic index and the differential line width (CRX and dLW, see Zechmeister et al. 2018). We also investigated the cross-correlation function (CCF, see Lafarga et al. 2020; Reiners et al. 2018a) to derive the full width at half maximum (FWHM), contrast (CON), and bisector span (BIS), and we derived a large number of additional indicators (see Schöfer et al. 2019). We searched for periodicities of all these indicators using the GLS periodogram. Many indicators show significant long-term signals around 365 d, and its 1 d aliases. The occurrence of this period in activity indicators of several other stars of our survey, in particular, of the other two targets discussed in this work, and the factthat it is compatible with one yearly cycle, makes it unlikely that stellar activity is the origin. This periodicity might be caused by small yearly environmental changes on the instrument or micro-tellurics that might affect the spectral line shapes to which the CCF and the measured pEWs are more sensitive than the actual RV measurements.

Because this yearly signal is not believed to be of stellar activity, and most importantly, because is far away from the planetary periods and derived stellar rotational period, we subtracted it so that we would be more sensitive to periods in the high-frequency regime. We show the residual GLS periodogram and its s-BGLS periodogram in Fig. 3. We found a signal at 121.0 d that is within the 1σ uncertainty of the photometric rotation period in TiO at 8430 Å with an FAP < 10⁻². Within the uncertainty of the first harmonic of the rotation period, we observed a peak in Hα with an FAP reaching almost 10⁻³. From the s-BGLS, the star showed the strongest activity in most indicators at periods attributed to the stellar rotation, whether at 120 d or 60 d, between January 2019 and October 2019 (CARMENES observation numbers 130–180).

Recently, signals at approximately 90 d in TiO 8430 Å and around 45 d in TiO 7050 Å have become significant. It is not clear where these signals originate. Recent works, for example, Shapiro et al. (2020) and Nava et al. (2020), have shown that the interplay of activity signals that is due to the distribution and different lifetimes of starspots and faculae on the stellar surface, may result in signals that cannot be directly attributed to the stellar rotation. However, we find a good agreement between our photometric results and spectroscopic results (see further down). The measured median pEW of the Hα line is + 0.00 ± 0.01 Å, and indicates that GJ 215 is not a Hα active star (Jeffers et al. 2018; Schöfer et al. 2019), which is in line with the long rotational period derived for this star, as is the upper limit of vsini < 2 km s⁻¹ measured by Reiners et al. (2018a).

Fig. 2 Joint GP model of the nightly binned photometric data of GJ 251. From top to bottom: T90 V, T90 R, TJO R, and SuperWASP.

5.3 Periodogram analysis and RV modeling

We show the results of the periodogram analysis of the CARMENES RV data in Fig. 4. A significant peak with an FAP < 10⁻⁷ is visible at 14.22 ± 0.01 d with an amplitude of 2.13 ± 0.23 m s⁻¹, as well as two additional peaks close to one day that we attributed to daily aliases of the 14 d period. Although the absolute GLS power and FAPs of the suspected aliases were smaller than the frequency of the 14.22 d signal, we used AliasFinder to verify that the 14.22 d signal represents the most probable true signal, which we confirmed. Additional strong secondary signals were identified at 73 d, and 119.5 d, each with an FAP of about 10⁻².

We also performed an independent periodogram analysis of the HIRES data. The strongest signal is at 604 d with an FAP of almost 10⁻⁵, followed by one at 14.2 d with FAP < 10⁻² (see Fig. 4). The latter is consistent with our strongest signal in the CARMENES data.

We combined the CARMENES and HIRES spectroscopic RV data by fitting an offset and jitter term for each instrument. The combined periodogram showed the highest peak at 14.24 d with an FAP < 10⁻¹¹, and the daily aliases of this signal were the second and third highest signals. Within our activity indicators, we did not identify any significant GLS periodogram peak with an FAP < 10⁻¹ at the frequency of the 14.2 d RV signal. Asignal at 14.731 ± 0.026 d in the Ca II IRT₃ line reaches almost 1% FAP, but is still larger than and can be well separated from the 14.24 d signal within the resolution of the GLS periodogram over the observed time baseline. We fit a Keplerian model to the signal at 14.24 d. Thelog-evidence of the different model fits applied to the data sets is given in Table 3.

The residual periodogram of the one-planet Keplerian fit on the HIRES and CARMENES combined data shows several remaining significant peaks at periods of 73.02 d (FAP < 10⁻⁵), 68.15 d (FAP < 10⁻⁴), and 67.86 d (FAP < 10⁻²). These periods are close to half of the derived rotational period of the star. We also observed a peak at 118.78 d with FAP < 10⁻², which is very close to the rotation period derived from photometry. As a simple test, we fit a sinusoid to the 118.78 d period. The signal at 67.86 d was then the most significant. It was necessary to fit an additional sinusoid for the 67.86 d signal to obtain a periodogram that did not show any signal with an FAP < 0.01, which showed that the other signals were connected through aliasing. The necessity of modeling two sinusoidal functions with periods close to the rotational period and its half suggests that these signals are caused by stellar activity, for instance, a multi-spot pattern, or amplitude variations caused by decreasing spot areas. To rule out the possibility of independent planet signals, we analyzed the coherence of these signals.

We used the s-BGLS periodogram to assess the coherence of the significant RV signals with increasing numbers of observations. We show the resulting s-BGLS diagrams in Fig. 5. We identified that neither the signals around 120 d nor the forest of signals between 50 d, and 73 d were stable over the observational time baseline. The 73 d signal lost about three orders of magnitude in signal probability after roughly observation 160 (May 2018), but reappeared in observation 210 (January 2019). All these mentioned signals showed a lack of coherence in the latest observations between observation 200 (January 2019) and 250 (October 2019). In contrast to these signals, the suspected planetary signal at 14.2 d never showed strong dips in its probability during the time of observations. These results together with the analysis of the activity indicators mean that this signal probably is of planetary origin. We refer to it as GJ 251 b.

Fig. 3 Generalized Lomb-Scargle periodograms of several activity indicators of GJ 251 from CARMENES spectroscopic data. The dashed black periodograms represent the GLS of the activity indicators, and the solid GLS periodogram represents the residuals after subtracting a 365 d signal. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The red dashed lines mark the rotation period and the first harmonic estimated from photometric data, while the dotted red lines show the 3σ uncertainties. The dashed black line marks a significant HIRES signal around 600 d, and the yearly period of 365 d. The dashed yellow line marks the period of the planetary signal published in this work.

Fig. 4 Generalized Lomb-Scargle periodograms of RV data for GJ 251 for CARMENES, HIRES, and a combination of both. The stellar rotational period derived by photometry is plotted as the dashed red line, and the 1σ and 3σ uncertainties are highlighted in red. We also indicate the harmonic of the rotational period and its uncertainty. The green line marks the suspected planetary signal. The blue lines mark the periods of the published planetary candidates by Butler et al. (2017) at 1.7 d, and the significant HIRES signal around 604 d.

5.4 A second planet in the system?

We found no evidence in the CARMENES and HIRES data for the planetary candidate claimed by Butler et al. (2017) at 1.74 d. Neither did we observe the 604-d signal in CARMENES data, which was highly significant in the HIRES data. We used the s-BGLS to verify whether these signals were more significant in the past but might have decayed over the time of observations. While we did not find any indication of the 1.74 d signal during the period of our observations, the 604 d signal showed variability in its signal probability (see Fig. 5). The 604 d signal already slightly lost coherence during the last HIRES observations. However, especially the CARMENES data show noncoherence of the signal. The fluctuation is a strong indication for a nonplanetary origin (Mortier & Collier Cameron 2017).

We performed a statistical test using model comparison in the framework of Bayesian evidence with juliet. We compared one-planet (P = 14.24 d) to two-planet models. Our results showed that the two-planet model with periods of 14.24 d and 1.74 d is not supported by the individual data sets or by their combination because its log-evidence is weaker than that of the simpler one-planet model. Fitting the 600 d signal as a second planet resulted in a significant model improvement compared to the one-planet model alone for the HIRES data. The same two-planet model (14.24 d, and 604 d) fit to the CARMENES data brought no significant improvement either compared to the one-planet model. Fitting the two-planet model to the combined CARMENES and HIRES data resulted in almost the same log-evidence as the one-planet model. Additionally, the derived planetary period at d deviates significantly from the d obtained from the fit on the HIRES data, even though we used a prior with an informative Gaussian distribution, hereafter referred to as normal prior, with mean of 604 d and σ = 30 d. These values were informed by the GLS periodogram peak in the HIRES data and its 3σ uncertainty.

We searched for any additional planetary signal hidden behind the stellar activity by once sampling a second Keplerian with a log-uniform prior between 15 and 8000 d and then sampling with a log-uniform prior between 0.5 and 14 d, while simultaneously modeling the stellar activity with a constrained GP model (see Sect. 5.5). We divided the two-planet model search into two runs for technical reasons: juliet needs a chronological order of the planetary periods. The two two-planet models combined with the GP showed no significant improvement compared to the one-planet and GP combined model, that is, the data suggest that all periodic variations except for the 14.2 d period are better or equally well described by a GP.

For the 600 d signal, we identified a periodicity with an FAP < 10⁻² in the Hα indicator at a period of 660 ± 21 d. The Na I doublet lines and the CRX showed a significant peak with an FAP < 10⁻³ at 300 d. The s-BGLS of the dLW shows that a signal close to 600 d was more significant in past observations around observation 150, which corresponds to April 2018 (see Fig. 3 again). Comparing the activity s-BGLS of the dLW to the s-BGLS of the RV data showed that this is about the same time at which the 600 d signal was most significant in the RV data. These results, along with our photometric results, suggest that the signals at 73, 119, and 600 d are not caused by Keplerian motion.

Fig. 5 Stacked-Bayesian GLS periodogram of the planetary signal at 14.22 d and on the zero-planet residuals (left diagram with its ownnormalization). The three s-BGLS on the right were calculated using the one-planet Keplerian RV residuals and show signals that we attributed to stellar activity: the forest of signals between 60 and 75 d corresponding to roughly half of the rotation period (left), the RV signals around the photometrically derived rotation period at 119.5 d (middle), and the long period signal around 600 d. In all four plots, the dashed black line indicates the boundary between HIRES and CARMENES data, which were taken successively.

Fig. 6 Posterior distribution of the GP fit to the RV data in the αGP vs. Prot plane for GJ 251. The color-coding shows the log-likelihood normalized to the highest value in the posterior sample. Gray samples indicate solutions with a Δ ln L > 10 compared to the best solution. Top: GP fit to the RV data with a wide uniform prior for the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric GP results and additional constraints on the other hyperparameters. We overplot the derived rotational period of the photometric GP and its 3σ uncertainties with vertical lines. The horizontal line marks the cut in αGP used to constrain the GP fit shown in the lower plot.

5.5 Simultaneous Keplerian and GP modeling

We performed a simultaneous fit of a one-planet Keplerian model together with the QP GP (Eq. (1)) to account for activity-induced RV variations. For the first GP model we used wide uninformative priors, which are shown in Table A.3, while we kept the same planetary and instrumental priors as for the one-planet fit (Table A.2). The posterior samples of this unconstrained GP can provide indications for the stellar rotational period given only the RV data (see also Angus et al. 2018; Stock et al. 2020).

Including the GP as a model for activity significantly improved the log-evidence (Δ logZ = 41 on the combined CARMENES + HIRES data) compared to the one-planet fit alone. We show the α_GP versus period diagram of the unconstrained GP posterior samples in the left plot of Fig. 6. Around 125 d we identified a region of higher posterior density, higher likelihood, and lower values of α compared to the rest of the posterior solutions. This indicates a more strongly correlated periodic signal. The derived median GP rotational period based on the CARMENES and HIRES RV data is d, which is consistent with the results from photometric data. For the final one-planet and GP simultaneous fit to the combined data, we applied several additional constrains on the GP to reduce the posterior volume of the model, which could lead to fitting incorrect residual signals (Angus et al. 2018). We applied a normal prior to the GP rotational parameter based on the stellar rotation period derived from the photometry and its 3σ uncertainty.

For Γ, which can be interpreted as the overall number of inflection points per function period, we applied a log-uniform prior between 10⁻¹ and 10¹. This prior is consistentwith about one to three local maxima per rotation period. Jeffers & Keller (2009) showed that this assumption is to first approximation valid for any stellar surface, independent of the number of starspots and their distribution. Similar, but even more informative priors on Γ, have been applied in several studies that used the QP GP kernel (see Nava et al. 2020, and references therein). The timescale parameter of the QP GP kernel, α_GP, is crucial formodeling a meaningful rotational signal. For instance, if it is large, then the squared exponential term of the kernel dominates, which allows for good fits to the data without requiring any periodic covariance structure, even when the data show clear periodicities (see also Angus et al. 2018). This effect is visible in the left plot of Fig. 6, where a plateau of posterior samples at high α_GP values populates the entire prior volume of the GP rotation parameter. As discussed in Sect. 5.1, it should not be strictly assumed that photometric and RV GP timescales are similar. A prior on α_GP based on photometry, as used for example for the stellar rotation, therefore needs future verification. For the moment, the upper boundary of the α_GP prior needs to be assessed for each target individually. Angus et al. (2018) proposed that this hyperparameter should be larger than the observed stellar rotation period. In the case of GJ 251, with a derived photometricrotation period of about 120 d, the suggested rule by Angus et al. (2018) would translate into α_GP < 3.5 × 10⁻⁵.

However, it is not clear whether this general rule can be applied to slowly rotating stars like GJ 251 because Angus et al. (2018) did not discuss such stars. For example, because of active longitudes (Jeffers & Keller 2009), a meaningful QP signal might still be detected that would be caused by starspots that decay over approximately half of the stellar rotation every time the active region points toward the observer. More importantly, from the unconstrained GP posterior samples, we find that a constraint on α_GP based on the rule by Angus et al. (2018) would mean that the overdensity of high-likelihood posterior samples detected at 120 d, given the data, would not be included in the final activity model. As expected, a GP model using the upper boundary of Angus et al. (2018) led to a log-evidence of − 710.2, which is about four lower than the unconstrained GP. Finding the right mixture between physical priors and data-driven behavior is critical for modeling stellar activity with GPs.

We constrained the upper boundary to α_GP < 3 × 10⁻⁴ d⁻². This constraint removed the plateau of posterior samples that fit noise on short timescales, which cannot be attributed directly to the stellar rotation and is captured by the instrument jitter parameter in our case. However, the observed high-likelihood posterior sample overdensity at the derived photometric rotation period is included in the GP model. We have applied similar constrains of α_GP with success in Stock et al. (2020). The priors of our constrained GP are summarized in Table A.3. The distribution of the posterior parameters of the constrained GP in the α_GP versus P_GP diagram is shown in the right plot of Fig. 6. Corner plots of all the fit parameters are provided in Figs. B.1 and B.2. The derived timescale parameter in our final GP model is d⁻², which translates into a decay time of d and is close to half the stellar rotation. The GP semiamplitude is consistent within 1σ with the expected RV semiamplitude due to stellar rotation estimated based on the relations of Suárez Mascareño et al. (2018).

In Fig. 7, we show our final median GP model. We calculated the GLS periodogram of the GP model to assess its temporal behavior. To calculate the GLS periodogram we chose to sample the GP model in two different ways: first, sampling identical to that of the original data, and second, sampled daily over the entire observation time. The GP model sampled as the real observations includes the true window function of the data. A visual inspection of this GLS periodogram shows that the highest peak is at 73 d, followed by another peak at 68 d. These were the most significant signals in the residuals of the simple one-planet fit. A peak at 28 d, about twice the planetary period and close to the lunar cycle, is also visible. The GLS periodogram of the GP model, sampled once a day over the entire observation time, shows that the GP does not model the 28 d period. The peak can be explained by the convolution of the GP model with the window function of the observations.

From the daily sampled periodogram decomposition of the GP models, the tuned GP does not model any periods close to the planet (or twice the planetary period). It models primarily the activity related signals at 63 d, 122 d, and 600 d. The unconstrained GP modeledthe 73 d signal more prominently than the signals at 120 d and 63 d, while the constrained GP modeled 120 d and 63 d more strongly. Nevertheless, the tuned GP was capable of producing the same strong peak at 73 d, given the data. This result shows in practice the conclusions and caveats given by Nava et al. (2020) that QP models can contain signals “unrelated to their true period”. The 73 d signal can be explained by an alias based on a sampling frequency of ~ 365⁻¹ d⁻¹ of the first harmonic at 63 d of the 120 d rotation period.

Finally, we show the combined Keplerian and tuned GP fit to the RV data, and a plot phased to the orbital period of GJ 251 b in Fig. 8. We display the final posterior solution of the planetary and GP parameters in Table 4. We derived a semiamplitude of m s⁻¹, a period of P = 14.238 ± 0.002 d, and an eccentricity of . The eccentricity of the system is consistent with zero because fits without this parameter provided similar log-evidence with fewer parameters. Based on our posterior samples and the stellar parameters (see Table 2), we derived further planetary parameters, which we also show in Table 4. According to this analysis, GJ 251 b has a minimum mass of 4.00 ± 0.40 M_⊕ and a semimajor axis of au.

Fig. 7 Gaussian process model for the RV data of GJ 251 without the planetary model, which is subtracted from the RV data. The constrained GP model is shown in red. The blue regions show 1σ, 2σ, and 3σ uncertainties. We show a zoom to some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP level of0.1%. In addition, we show the unconstrained GP model as the dashed black line in the upper plots and as the gray periodograms in the lower plots.

Fig. 8 Left: radial velocity data with a combined model of one-planet and stellar activity using a Keplerian model and a GP. Right: plot phased to the orbital period of GJ 251 b without a GP component. The bottom plots show the residuals after the fit.

5.6 Transit search and analysis with TESS

GJ 251 was observed with the TESS satellite (Ricker et al. 2015) in sector 20, in the period from 24 December 2019 to 21 January 2020, with a total of 16 556 data points, but was not marked as a TESS object of interest (TOI). We independently searched for a transit signal using the transit least-squares (TLS; Hippke & Heller 2019) algorithm on PDCSAP light curve. We did not identify any TLS signal that could be attributed to any possible transit for GJ 251. However, we identified a significant sinusoidal-like signal with a frequency f ≈6 d⁻¹ (period 0.165 d) and (1.96 ± 0.4) × 10⁻⁴ relative fluxamplitude variation, equivalent to about 0.2 mmag), as well as its harmonics f∕2 and 2f with lower amplitudes. We show a phase plot of this signal in Fig. 9. This signal is not observed within the RV data of this target.

This 4 hsignal was already present in the simple aperture photometry (SAP) light curve. We checked that it was not instrumental in origin by extracting and analyzing the light curves of the 968 objects present in the same TESS S20 sector, Camera 1, and CCD 3 as GJ 251. No other star showed the same periodicity.

We considered the possibility that the signal originated from thermodynamical excitations of p- and g-modes, as theoretically predicted by Rodríguez-López et al. (2014). However, the pulsation hypothesis is unable to explain the presence of subharmonics of the main frequency in the periodogram. Moreover, solar-like pulsations or granulation, which have not yet been detected in M-dwarf stars, can also be discarded because they are predicted to be on the order of minutes.

Lucky-imaging observations with FastCam (Cortés-Contreras et al. 2017) and Robo AO images (Lamman et al. 2020) have not detected any resolved visual companion. The TESS aperture includes several objects. In particular, the two brightest stars in the aperture mask are only 4.87 mag and 5.95 mag fainter in the G band, corresponding to a flux contribution of 1.1% and 0.4%, respectively. A 18 mmag sinusoidal amplitude variation in the former or a 50 mmag amplitude variation in the latter could account for the detected 4 h signal. We chose different subapertures to extract the light curve from different regions of the TESS full-frame images. It did not affect the amplitude of the short-period signal, making it unlikely that the periodicity originated in background contamination. With our analysis, we cannot draw any final conclusion on the origin of the 4 h signal for GJ 251.

We estimated the radius of GJ 251 b with the mass-radius relation of Zeng et al. (2016) and assumed an Earth-like core-mass fraction of 0.26 to be approximately 1.48 R_⊕, which translates into a transit depth of roughly 1.4 ppt for GJ 251 b. Such a signal should be detectable by TESS in case of a full transit. However, we were unable to detect any transit in the light curve, in particular given the estimated t₀ and the orbital period P of GJ 251 b by the RV fit and their uncertainties. In particular, we ruled out a transit event within 1σ of t₀, but not within 3σ, because of an observational gap in the TESS light curve. Unfortunately, GJ 251 will not be observed again by TESS.

6 HD 238090

6.1 Photometric monitoring

We took ground-based photometry for HD 238090. We combined our T90 and TJO data with public data from MEarth taken in 2009, 2010, and 2014, as well as data from NSVS taken between May 1999 and March 2000. A periodogram analysis of the combined data sets indicated periodicities around 100 d. Using juliet, we fit for an offset and jitter terms for each instrument and filter. We used the same GP kernel and priors as for the photometric analysis of GJ 251 and separated the amplitude hyperparameters for each instrument while keeping global hyperparameters for the timescale and rotation. The priors are given in Table A.1. From our GP analysis, we derived that HD 238090 has a rotation period of d. Figure 10 shows the distribution of the posterior samples in the α_GP -P_rot. space, and Fig. 11 shows the median GP model of each photometric data set together with the data and uncertainties. Following the same approach as for GJ 251 and applying the relations by Suárez Mascareño et al. (2018), we derived the median and expected semiamplitude to (mean at − 5.69) and (mean at 4.25 m s⁻¹).

Fig. 9 TESS PDCSAP light curve of GJ 251 folded to the 0.165 d signal. The black dots and uncertainty bars represent binned TESS data.

6.2 Spectroscopic activity indicators

The periodogram analysis of our spectroscopic activity indicators is displayed in Fig. 12. We identified signals with an FAP < 10⁻³ in some indices close to 1 d and 365 d (see Sect. 5.2 for a discussion of these signals). After subtracting the yearly signal, the GLS periodograms of the residuals of many indicators show a signal around 480 d with FAP < 10⁻³ and a long-term trend. We also identified a signal in the FWHM CCF with an FAP < 10⁻² at 106.1 d, which is close to the derived stellar rotation period. Within the s-BGLS of the residual activity indicators of TiO, which we show in Fig. 12, we found various signals around 50 d, which is roughly half the photometrically derived stellar rotation period. These signals were more significant in previous observations between July 2018 and February 2019 (CARMENES observations60 to 90). Overall, the activity indicators show that the star exhibits no significant level of activity at periods shorter than 20 d over the time of RV observations.

Fig. 10 Posterior distribution of the GP fit to the photometric data in the αGP vs. Prot plane for HD 238090. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10 compared to the best solution.

6.3 Periodogram analysis

The GLS periodogram of the RV data for HD 238090 is shown in Fig. 13. A significant peak with an FAP < 10⁻⁵ is visible at 13.68 d. Two additional signals with almost the same GLS power accompany this signal at periods close to 0.93 d and 1.08 d. No additional signals were significant in the data. The investigation with AliasFinder led to the conclusion that the 13.68 d period is the most probable true period of the sampled signal because simulated periodograms based on this period fit the observed periodogram better than the periods close to one day. We show the corresponding plot obtained by AliasFinder in Fig. C.1.

An investigation of the 13.68 d signal with the s-BGLS showed that the signal was coherent and increased in probability over the observation time. We display the s-BGLS of the signal in Fig. 14.

6.4 RV modeling

Becausethe signal at 13.68 d showed long-term coherence and had no counterparts in the activity indicators, we fit a Keplerian model to the 13.68 d signal, hereafter referred to as HD 238090 b. Table 5 shows that a one-planet fit is significantly favored by the data compared to a flat model.

The residual periodogram of a one-planet fit to the 13.68 d signal showed peaks with an FAP of almost 10⁻² for 106.4 d. The signal at 106.4 d is below our optimal detection criterion by the GLS analysis, which is FAP < 10⁻³. Additionally, it resides within the 3σ uncertainty of the derived rotation period of the star by photometry and has a counterpart in the FWHM CCF, as discussed before. An s-BGLS analysis of the 106 d signal indicated that the probability of the signal decreased by almost two magnitudes after 30 observations before reappearing in later CARMENES epochs. We show the s-BGLS in Fig. 14.

These results indicate that the signal is probably of non-Keplerian origin. Nevertheless, we performed an additional statistical model comparison where we compared a second circular Keplerian to a GP for this signal. The two-planet fit resulted in a log-evidence improvement compared to a one-planet model of Δ ln Z = 5.2. A wide unconstrained GP to account for the 106 d signal performed equally well. Statistically, there is no clear tendency for a GP or Keplerian model. Future observations of HD 238090 might be warranted to completely exclude the possibility of a Keplerian signal with a period of 106 d, especially because such a planet would reside inside the optimal habitable zone described by Kopparapu et al. (2013). Based on the current data and our photometry, s-BGLS, and activity indicator analysis, we regard the 106 d signal as caused by stellar activity.

As for GJ 251, we improved the GP modeling by constraining the prior volume (and therefore the posterior). The GP alpha-period diagram (α_GP versus P_GP) of the unconstrained GP is shown in the top plot of Fig. 15 and the priors are given in Table A.3. We found a posterior overdensity with a marginally higher likelihood around 50 d, which is close to half the derived photometric rotation period for the unconstrained GP. However, the distribution of posterior samples showed no peculiar structure overall and suggested that the RV data of HD 238090 are not significantly affected by a strong correlated quasi-periodic signal with decay-timescales of more than several days.

In the lower plot of Fig. 15, we show a GP for which we constrained the rotation parameter P_GP to be Gaussian distributed around the derived stellar rotation period and the timescale parameter α_GP to be between 10⁻³ and 10⁻²⁰. With the new priors, the dynamic nested sampling algorithm found posterior solutions around 106 d that reached a similar maximum likelihood as the unconstrained GP posterior samples. The log-evidence of this model was equal to the unconstrained GP. That the distribution of α_GP reached the prior boundary at 10⁻²⁰ implies that this parameter converged to zero because this decay timescale is orders of magnitudes longer than the observation time.

These GP results motivated us to fix the α value to 10⁻²⁰, which is consistent with zero. This resulted in only three free hyperparameters to be fit by the GP. This simpler rotational GP performed similarly in terms of log-evidence as the two previously described QP-GPs. We applied this GP as the final activity model of the system, hereafter called constrained GP, and show the model in Fig. 16. For comparison, we plot the unconstrained quasi-periodic GP model, which included high α_GP values within its posterior, in the same figure. The constrained GP model represents a more realistic fit of the data. The GP semiamplitude is consistent within 1σ with the expected semiamplitude that is due to the stellar rotation, estimated based on the relations of Suárez Mascareño et al. (2018), and the GP rotation period is within the 3σ uncertainty of the photometric estimate of the stellar rotation. The GLS periodogram decomposition of the constrained GP model sampled daily indicates that the GP only models the 106 d period and its first harmonic.

We performed a search for any additional planetary signal hidden behind the stellar activity with the tuned activity model. For this, we included a second Keplerian signal using a log-uniform prior between 0.5 d, and 13 d and then 14 d to 1000 d for its period, while simultaneously fitting for GJ 251 b and the stellar activity with the constrained GP. We found that the two two-planet models together with the GP did not perform significantly better than the one-planet plus GP model. We therefore preferred the one-planet model for its simplicity.

We show the final one-planet and activity model to the RV data in Fig. 17 and the posterior parameters in Table 4. We derived a nonsignificant eccentricity of , as a circular model resulted in similar log-evidence. We derive the minimum mass of HD 238090 b to be M_⊕. We provide corner plots of all posterior samples in Figs. B.3 and B.4, respectively.

Fig. 11 Gaussian-process model for each nightly binned photometric data set of HD 238090. From top to bottom: NSVS, MEarth 2009 to 2010, MEarth 2010 to 2014, SNO V, SNO R, and TJO. For each instrument, we fit individual GP hyperparameters for the amplitudes and Γi, but we used global GP hyperparameters for the timescale of the amplitude modulation and the rotation period.

Fig. 12 Generalized Lomb-Scargle periodograms of several activity indicators based on spectroscopic data obtained by CARMENES for HD 238090. The dashed black periodograms represent the GLS on the activity indicators, and the solid GLS periodogram represents the residuals from which a 365 d sinusoidal signal was subtracted. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The solid red lines mark the rotation period estimated by photometric data, and the dashed red lines show the 3σ uncertainties. The dashed black line marks a period of 365 d, and the dashed yellow line marks the period of the planetary signal published in this work.

Fig. 13 Generalized Lomb-Scargle periodogram of RV data of HD 238090 of the zero-planet fit, one-planet fit, and one-planet GP simultaneous fit. The stellar rotational period derived by photometry is plotted as the dashed red line, and 1σ and 3σ uncertainties are highlighted in red. The green line indicates the planet period at about 13.7 d.

Fig. 14 Stacked-Bayesian GLS of the planetary signal (left) and a signal close to the estimated stellar rotation (right).

Fig. 15 Posterior distribution of the GP fit to the RV data in the αGP vs. Prot plane for HD 238090. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10. Top: GP fit to the RV data with a wide uniform prior to the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric GP results and an upper αGP constraint.

Fig. 16 Gaussian-process model for the RV data of HD 238090. The planetary signal is not included in the model and subtracted from the RV data. The constrained GP model is shown in red. The blue regions shows the 1σ, 2σ, and 3σ uncertainties. We show a zoom into some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP of 0.001. We also show the unconstrained GP model as the dashed black line in the upper plots and the gray periodograms in the lower plots.

6.5 Transit search with TESS

TESS observed HD 238090 (TIC 224289449) in sectors 15, 21, and 22. Based on the mass-radius relation by Zeng et al. (2016) and applying a core-mass fraction of 0.26, which corresponds to Earth-like composition, we estimated a planetary radius of 1.69 R_⊕ for HD 238090 b. Given the stellar parameters of HD 238090, the transit depth was approximated to 0.72 ppt, which should be detectable in the TESS light curve. We investigated the TESS light curve around the estimated t₀ from the RV fit. We did not identify any transit event.

We calculated the TLS periodogram of the light curve and found three signals with a signal detection efficiency (SDE) > 7, which we regard as significant (Hippke & Heller 2019), at periods of 30.00 d, 27.38 d, and 13.67 d. The first periodsare close to the TESS sector length, while the 13.67 d period is about half the TESS sector length, but represents exactly the period of HD 238090 b. However, the time of transit center derived for the 13.67 d signal detected in the TLS is incompatible with the value obtained from the RV analysis (t_0,transit = 2458638.30 ± 0.05 BJD versus t_0,RV = 2458630.35 ± 0.07 BJD). Inspection of the TESS light curve showed that the presumed transits were fit at the edges or inside observational gaps of the light curves of the TESS sectors. We concluded that these TLS signals, although significant and close to the planetary period, represent false positives caused by the observational sampling and that no transits of HD 238090 b are detected.

7 Lalande 21185

Díaz et al. (2019) reported the discovery of a temperate super-Earth orbiting Lalande 21185 with a period of 12.95 d or, less likely, 1.08 d, but the period could not be determined unambiguously because of aliasing. Díaz et al. (2019) did not find evidence for a planetary candidate reported previously by Butler et al. (2017), which was supposed to orbit the star at 9.9 d. We reanalyzed the HIRES and SOPHIE data together with our CARMENES observations.

7.1 Photometry and spectroscopic activity indicators

Díaz et al. (2019) took extensive photometric observations between 2011 and 2018 with the Tennessee State University T3 0.40 m automatic photoelectric telescope at Fairborn Observatory in southern Arizona. Based on their analysis, they reported a stellar rotation period of 56.15 ± 0.27 d for Lalande 21185. This rotation period is consistent with values previously obtained by Noyes et al. (1984, 48 d) and by Oláh et al. (2016, 54 d), although these studies did not provide uncertainties.

Based on the derived stellar rotation period by Díaz et al. (2019), we made use of the relations by Suárez Mascareño et al. (2018) and derived for the median and median expected RV semiamplitude (mean at − 5.37) and (mean at 7.27 m s⁻¹), respectively.

The periodograms of the activity indicators are provided in Fig. 18. After subtracting the yearly signal (see Sect. 5.2 for discussion of this signal), we found several signals in the CARMENES activity indicators close to the reported stellar rotational period of Lalande 22185, e.g., TiO at 7050 Å with an FAP < 10⁻¹, TiO at 8430 Å and at 8860 Å with an FAP < 10⁻², in the FWHM with an FAP < 10⁻¹ and in BIS with an FAP < 10⁻². The s-BGLS of these signals show the instability of these signals over the observation time. For the activity analysis of the SOPHIE data, we refer to Díaz et al. (2019).

In addition to these stellar rotation related signals, we observed two significant signals (FAP < 10⁻³) in Hα at periods of 1313 d and 539 d, and around 1400 d in the dLW, the Na lines, and the contrast CCF. Additionally, we observed a linear trend in the CRX. In the CRX, dLW, and the Na lines, we identified signals at roughly 14 d (with FAPs of <0.1, <0.01, and ~0.1, respectively), which is about 1∕4 of the stellar rotational period, but close to the claimed planetary signal by Díaz et al. (2019). Owing to the long CARMENES time baseline, the GLS periodogram resolution allows us to separate these peaks from the planetary signal at 12.95 d.

Fig. 17 Left: radial velocity data with a combined model for one planet and activity using a Keplerian model and a GP. Right: plot phased to the orbital period of HD 238090 b without a GP component.

7.2 Periodogram analysis

Figure 19 displays our GLS periodograms for the CARMENES, SOPHIE, and HIRES data, as well as residual periodograms for different models fit to the combination of CARMENES and SOPHIE data. The CARMENES and the SOPHIE GLS periodograms show a significant peak at a period of 12.95 d, which is the period of the planetary signal published by Díaz et al. (2019). With our CARMENES data alone, we can confirm this signal. In addition to the 12.95 d signal, the SOPHIE data show a signal at 55 d, with an FAP < 10⁻², which is consistent with the stellar rotation period of 56.15 ± 0.27 d. The CARMENES data also show a significant long-period signal at 2677 d. A linear trend has been reported by Díaz et al. (2019) from the SOPHIE data. With the addition of the CARMENES data, we now observed a long period.

Neither the CARMENES data nor SOPHIE data or their combination shows a significant signal at 9.9 d, where a signal was claimed by Butler et al. (2017) using HIRES data. This signal is visible with an FAP < 10⁻³ in the periodogram of the HIRES data, but we find many more signals with similar or higher significance in the periodogram of the HIRES data, with amplitudes of a few m s⁻¹. All these signals are absent in the SOPHIE and CARMENES data, however. The time baseline of the combined CARMENES and SOPHIE observations is about 8.2 yr, and the precision of both instruments should be appropriate to identify these signals if they were still present in the RVs of Lalande 21185. The highly significant presumably planetary signal at 12.95 d, with an amplitude of roughly 1.4 m s⁻¹ in the CARMENES and SOPHIE data, cannot be identified in the periodogram of the HIRES data. The amplitude of the signal might not be large enough for HIRES because it is at the limit of the long-term precision, which has been about 1–2 m s⁻¹ since 2004 (Butler et al. 2017). The sampling of the HIRES data shows many nights with multiple observations. After applying a nightly binning scheme on the HIRES data, the GLS periodogram showed no peak with an FAP < 10⁻². Because the 12.95 d signal is absent in the HIRES data and the forest of significant but spurious signals at various frequencies, we restrict the RV analysis to the combined SOPHIE and CARMENES data sets and treat the HIRES data individually.

Díaz et al. (2019) reported that the SOPHIE data did not allow determining the orbital period of the planetary signal unambiguously because of aliasing. By adding our CARMENES data and using the AliasFinder, we can confirm that the sampled signal has a period of 12.95 d because the AliasFinder simulations were able to reproduce the properties of the observed periodogram only when this period was assumed to be the correct one. We show the relevant plots in Fig. 20.

7.3 RV modeling with CARMENES and SOPHIE

We fit a Keplerian model to the 12.95 d signal (see Table 6 for the log-evidence). The GLS periodogram of the residuals of the combined data set reveals additional significant signals at periods of 2852 ± 568 d with an FAP < 10⁻³ and at 55.3 ± 0.2, 61.3 ± 0.3, and 383.8 ± 10.6 d with an FAP < 0.01. The long-period signal is highly significant in the combined data set, but we find significant variations at similar periods or half this period in some of the activity indicators, mainly in CRX and Hα.

The s-BGLS was used to assess the coherence of all these signals. We show the s-BGLS of the combined data set around the orbital period of the planet signal at 12.95 d and the s-BGLS of the one-planet fit residuals around the long-periodic signal and the stellar rotational period in Fig. 21. The data are ordered chronologically. The signal probability of the suspected planetary signal at 12.95 d increases and shows coherence over the entire observation time. The signal at 55 d, which we related to the stellar rotational period, shows an increase in signal probability until roughly 360 observations and has decreased since then by about four orders of magnitude in probability. A similar pattern but anticorrelated to the rotational signal in terms of signal probability over time is visible for the long-period signal.

Based on the results on the activity indicators and the s-BGLS analysis, we find that the period around 2800 d can be best explained by a long-term activity cycle. We fit a simple sinusoid to this signal in order to search for additional signals. The residual periodogram of the one-planet + sinusoid fit has one remaining significant signal at 55.3 d. A fit of this signal, which represents the rotational period, with a second sinusoid, resulted in a flat periodogram; no peaks in the GLS periodogram reach an FAP < 0.1. Based on CARMENES and SOPHIE RV data, the Lalande 21185 system can be explained by one Keplerian model with an orbital period of 12.95 d, and two sinusoids that model the activity contribution.

We fit the system using a one-planet model simultaneously with a GP model, which accounts for these activity-related signals. In a first step, we fit a rather unconstrained GP to the data simultaneously with the one-planet Keplerian model. In the top plot of Fig. 22 we show the posterior sample distribution in the α_GP versus P_rot plane for this fit to Lalande 21185. Most of the posterior samples peak with a bimodal distribution at a period of 55 d and 65 d and at rather low α-values representinga stable quasi-periodic signal. We also see fewer posterior samples around 100 d but with a lower likelihood. The 65 d signal was already visible in the residuals of a simple one-planet fit and belonged to an alias of the rotational period based on a yearly sampling frequency f_s = 1∕365.25 d⁻¹. In the next step, we fit a more constrained GP with a normal prior on the GP rotational period based on the photometric estimates and its 3σ uncertainty and additional constraints in Γ and α deduced from the posterior distribution of the unconstrained fit as before. The priors of the applied GP models are provided in Table A.3. This GP fit resulted in a significant improvement of the log-evidence compared to the unconstrained GP and represented the best model we derived for this system (see Table 6). This GP fit is compatible with the estimated RV semiamplitude from the stellar rotation period given the relations by Suárez Mascareño et al. (2018).

We show the GP model and a GLS periodogram analysis to assess its temporal behavior in Fig. 23. We identify that in addition to the 56 d and long-term period, signals at about 380 d and 60 d are modeled. As for the other two targets, we searched for planetary signals hidden behind the stellar activity by sampling for a second Keplerian, using a log-uniform prior between 0.5 d, and 12 d and then a log-uniform prior from 13 d to 3000 d for the second Keplerian, while simultaneously fitting the stellar activity with our final GP model. The two two-planet models combined with the GP performed worse in log-evidence (see Table 6) than the one-planet model combined with the GP; this means that no additional Keplerian signal is statistically supported by the data. A plot of the final one-planet and GP model, the RV data, and a phase plot to the planetary period of Lalande 21185 b are provided in Fig. 24. The updated orbital parameters of Lalande 21185 b based on the posterior solutions of the fit to the combined CARMENES and SOPHIE data are given in Table 4. Additionally, we list further derived planetary parameters, such as the planetary minimum mass, which is estimated to 2.69 ± 0.25 M_⊕.

Fig. 18 Generalized Lomb-Scargle periodograms of several activity indicators based on spectroscopic data obtained by CARMENES for Lalande 21185. The thinner dashed black periodograms represent the GLS on the activity indicators, and the solid GLS periodogram represents the residuals from which a 365 d sinusoidal signal was subtracted. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The dashed red lines mark the rotation period and the first harmonic estimated by photometric data, and the dotted red lines show the 3σ uncertainties. The dashed black lines mark a period of 365, 1400, and 2800 d, and the dashed yellow line marks the period of the planetary signal published in this work.

Fig. 19 Generalized Lomb-Scargle periodograms of Lalande 21185 for the different data sets and the combined CARMENES and SOPHIE data. Residual periodogram of the one-planet fit and the one-planet GP simultaneous fit for the combined CARMENES and SOPHIE data are also shown. Marked frequencies represent the claimed planet at 12.95 d by (Díaz et al. 2019; green), photometric stellar rotational period (red), the claimed planet candidate by (Butler et al. 2017; orange), and the long-term period (blue).

7.4 HIRES RV data

We excluded the HIRES data from our final analysis of Lalande 21185 because of spurious frequencies, noisy data and the absence of an obvious planetary signal. We fit a one-planet model with the same priors as for the SOPHIE and CARMENES data to the HIRES data. For example, for the period, we used . While this improved the log-evidence significantly, it resulted in an inconsistent orbital period of d and an extremely high eccentricity of .

We also tested more informed priors. For instance, we applied normal priors to every planetary parameter based on the solution of the fit to the CARMENES and SOPHIE data. This fit performed reasonably well because it was still significantly better than a flat model, which could indicate that the planetary signal is apparent in the HIRES data. However, the same model fit to the daily-binned HIRES data did not result in a log-evidence improvement compared to a zero-planet model. The noise level of the HIRES data compared to the other data sets and the small planetary amplitude, which is at the limit of the HIRES long-term precision, justifies the exclusion of this data set. Nevertheless, we used the extended HIRES time baseline to analyze the long-period signal. With HIRES, we have a total of 737 RV observations for Lalande 21185. We find that the long-period signal is also apparent in the HIRES data. However, the HIRES data between BJD 2452200 and 2453200 indicate a possible phase shift of that signal, consistent with our analysis that this signal is caused by a long-term activity cycle. We fit our simplistic model of one planet at 12.95 d, and two sinusoids, to the combined CARMENES, SOPHIE, and HIRES data. The residual GLS periodogram of this fit had no peaks with an FAP < 10⁻³, indicating that the HIRES data do not indicate an additional coherent signal when combined with CARMENES and SOPHIE data.

7.5 Transit search with TESS

TESS observed Lalande 21185 in sector 22, but the star was not announced as a TOI. We ran our independentsignal search with the TLS on the PDCSAP light curve. We find only one peak with SDE > 7 (Hippke & Heller 2019), which is at 13.0 d, very close to the planetary orbital period derived from the RV fit, just as for HD 238090.However, peaks in this period range are visible in many light curves and are often caused by the observational gap of the TESS observations. Nevertheless, as for HD 238090, we performed a transit fit to this signal using juliet to evaluate how such a fit performs compared to a nontransit model. The log-evidence of a transit model was not significantly different (Δ logZ ≈ 0.5), indicating that a simpler nontransit model is better than a transit model. The derived minimum mass of the planet was used to estimate the radius with the mass-radius relation of Zeng et al. (2016). With an approximated radius of 1.33 R_⊕, the transit depth was approximated to 0.97 ppt for Lalande 21185 b. We were unable to identify any transit events around the RV estimated time of transit center for Lalande 21185 b.

Fig. 20 Alias tests for Lalande 21185. The top plot (a) shows simulations motivated by a sampling frequency of . The bottom plot (b) shows simulations motivated by a sampling frequency of fs2 = 1.0027 d−1. Each row in these plots corresponds to one set of simulations for which the frequency of the injected signal is indicated by a vertical dashed blue line. The first row shows simulations with a period of 12.95 d, and the second and third row show the simulations in which the first-order aliases of 12.95 d, regarding the investigated sampling frequency, were injected. Each column shows informative ranges of the periodograms based on the assumed sampling frequency and can be used to compare data and simulations. From 1000 simulated data sets, we show the median of the obtained periodograms (solid black line), the interquartile range, and the ranges of 90% and 99% (gray shades). The periodogram of the observed data is plotted with a solid red line. The angular mean of the phase and the standard deviation is shown in the clock diagrams (black line and gray shades) and can be compared to the phase ofthe signals in the observed periodogram (red line).

8 Discussion

8.1 GJ 251

Based on our analysis of CARMENES and HIRES RV data, we report the discovery of GJ 251 b, a planet that orbits its host star with a period of 14.24 d. The posterior sample median eccentricity is , which is notsignificant and consistent with zero, as shown by a log-evidence comparison with circular Keplerian models. In Fig. 25, we place GJ 251 into context with other confirmed exoplanets around M dwarfs detected with the RV method and the other two planets discussed in this work. The equilibrium temperature of GJ 251 b, assuming a zero Bond albedo, is K (see Table 4). Our minimum mass and temperature estimates add GJ 251 b to the family of temperate super-Earths. The planet receives about 2.5 times the flux of Earth. Of the three planets discussed in this work, GJ 251 b is the planet with the most moderate temperature. However, according to Kopparapu et al. (2013), the planet is too close to its host star to be in the habitable zone, that is, to allow liquid water on its surface.

Butler et al. (2017) previously claimed a planetary candidate for GJ 251 with an orbital period of 1.74 d. Our analysis of the CARMENES data and the combined data from CARMENES and HIRES did not confirm this claim. An independent analysis of the HIRES data does not identify any signal close to the 1.74 d period as significant. We note that Butler et al. (2017) used a different approach based on autocorrelation functions and a statistical model for the RVs, including a moving average to model correlated noise and information provided by the activity index based on Ca II lines.

The HIRES RV data show a significant second RV signal around 600 d. This signal was also visible during some epochs in the CARMENES data, but because the signal is incoherent and a significant Hα activity signal lies at the same period, we attribute this signal to nonplanetary origin. Any other additional significant signals can be best explained by the stellar rotation of GJ 251 or its harmonics. In the current RV data, we find no evidence for a second companion in the GJ 251 system. There is also no significant linear trend in the RV data over the entire time baseline of the combined HIRES and CARMENES observations, which is about 15.0 yr. However, the Gaia DR2 (Gaia Collaboration 2018) catalog lists a significant astrometric excess noise for this star, which could be caused by a massive companion on a wide orbit.

Fig. 21 Stacked-Bayesian GLS periodogram of Lalande 21185. The color bar to the left color-codes each observation with the associated spectrograph (orange: CARMENES, and teal: SOPHIE). Left: s-BGLS of the zero-planet model at the period of theplanetary signal. Middle: s-BGLS on the residuals of the one-planet fit with a period of 12.95 d shown around a period of 55 d, which is the rotational period determined by Díaz et al. (2019). Right: s-BGLS on the residuals of the one-planet fit around the observed long-period signal.

8.2 HD 238090

Our analysis of CARMENES RV data shows that HD 238090 is orbited by a warm super-Earth, HD 238090 b, with an orbital period of 13.69 d. The posterior sample median eccentricity is , which, as inthe case of GJ 251, is not significant. HD 238090 b has a minimum mass of roughly 6.8 M_⊕, and orbits its host star at a separation of approximately 0.093 au with an equilibrium temperature of 470 K (see Table 4). With total insolation about eight times that of Earth, the planet istoo close to the host star to sustain liquid water on its surface. Figure 25 shows the position of the planet in the minimum mass-period plane.

8.3 Lalande 21185

Our analysis of the CARMENES and combined CARMENES and SOPHIE RV data for Lalande 21185 confirms the findings by Díaz et al. (2019) regarding Lalande 21185 b. With our data, we can break the degeneracy between the daily aliases and confirm thatthe planet orbits the star with a period of 12.95 d. With the additional CARMENES observations, we have reduced the uncertainties of most planetary parameters by a factor of two or more compared to the estimates by Díaz et al. (2019).

The CARMENES data of Lalande 21185 agree well with the SOPHIE data but cast doubt on the HIRES data for this system. With 476 precise high-precision RVs for the system by CARMENES and SOPHIE, we find no evidence for a second planet in the system. In particular, the planet candidate claimed by Butler et al. (2017) at a period of 9.9 d is absent in the CARMENES and SOPHIE RV data. In addition to the 12.95 d signal, there are two additional significant signals in the combined CARMENES and SOPHIE data set at 55.3 d and 2800 d, which can be attributed to the stellar rotation and possible long-term activity, respectively.

Fig. 22 Posterior distribution in the GPα vs. Prot plane for Lalande 21185. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10. Top: GP fit to the RV data with a wide uniform prior to the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric rotational period proposed by Díaz et al. (2019).

Fig. 23 Gaussian process model for the RV data of Lalande 21185. The planetary signal is not included in the model and subtracted from the RV data. The constrained GP model is shown in red. The blue regions shows the 1σ, 2σ, and 3σ uncertainties. We show a zoom into some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP of 10−3. We also show the unconstrained GP model as the dashed black line in the upper plots and as the gray periodograms in the lower plots.

8.4 Formation scenario

The planetary systems presented in this work represent systems that share similar orbital properties (period, separation, insolation, and equilibrium temperature) and minimum masses that place them into the group of temperate or warm super-Earths. While the available data are insufficient to draw definite conclusions on the exact formation channel of these systems, the derived orbital and planetary parameters allow for some cautious conjectures. These planets are commonly thought to form by combined accretion of planetesimals and pebbles (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch 2019).

Their final masses can be best explained if the supply of pebbles was cut off during their formation, preventing their evolution into gas giants. One way to stop the supply of solid material from the outer disk is the emergence of a massive companion that stops the pebble flux by opening a gap in the protoplanetary disk (Ormel 2017). However, we do not see evidence for additional planets in any of the three systems, even though the available long-baseline data allow for a strong sensitivity for the detection of such companions in the cases of GJ 251 and Lalande 21185.

Another proposed mechanism to terminate pebble accretion is self-isolation from the pebble flux by modulation of the gas pressure profile by the growing planet (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014). For all three planets, their combination of mass and orbital separation would be consistent with this possibility. In this case, the measured masses are simply their pebble isolation masses, regardless of other planets in the system. On the other hand, systems such as those presented here can also be explained by models that grow solid cores only by accretion of planetesimals. In this scenario, planetary core growth can reach a natural stall before entering runaway gas accretion because the accretion efficiency is lower than pebble accretion (Emsenhuber et al. 2020a,b). In a related study focusing on planet formation around low-mass stars, planets similar to our discoveries are among the most abundant in a synthetic planet population (Schlecker et al. 2020; Burn et al., in prep.). While the model by Burn et al. (in prep.) predicts an average multiplicity higher than one for super-Earths, the singular detection in GJ 251, HD 238090, and Lalande 21185 could be explained by planets with lower masses or wider orbits that do not reach the RV semiamplitudes necessary for robust detection, especially in the case of strong activity.

Fig. 24 Left: radial velocity data with a combined model for one planet and activity using a Keplerian model and a GP. Right: plot phased to the orbital period of Lalande 21185 b without a GP component.

Fig. 25 Confirmed exoplanets around M-dwarf host stars, detected with the RV method and with minimum masses below 10 M⊕ and orbital periods between 1 d and 400 d. The color-coding represents the spectral type. GJ 251 b, HD 238090 b, and Lalande 21185 b are marked in red. Data taken from http://exoplanet.eu

9 Summary

We presented the discovery of two super-Earth planets around the low-mass stars GJ 251 and HD 238090 with orbital periods of 14.24 d, and 13.67 d, respectively, based on CARMENES VIS RV observations. For GJ 251, we additionally used RV data obtained by HIRES in order to increase the time baseline and to search for additional signals. We also confirmed the nearby temperate super-Earth Lalande 21185 b recently discovered by Díaz et al. (2019), and we can robustly determine its orbital period to be 12.95 d. No transits could be detected with TESS for any of the three systems. The RV data of GJ 251 and Lalande 21185 exhibit long-term periods, which we attribute to activity. Furthermore, all three systems show RV signals related to the stellar rotation period in the RV residuals of the planetary fits.

We modeled the stellar activity using GP models based on a quasi-periodic kernel simultaneously with the Keplerian signals. In particular, we carefully modeled the stellar activity by applying physically motivated constraints to the GP hyperparameters to ensure that the GP did not fit any signal unrelated to stellar activity. We advocate the use of classical periodograms to decompose the modeled frequencies by GPs. Such an analysis can be used as verification of the desired GP behavior on the data set. Nevertheless, the unconstrained GP posterior distribution can provide useful information on the stellar activity. In particular, we used the GP timescale versus GP rotation plane to infer more information about possible rotation periods and lifetimes of stellar surface features. For the analysis of the RV data, we used various advanced tools, such as juliet together with a Bayesian approach based on the log-evidence, AliasFinder to further distinguish samples signals and their aliases, and the s-BGLS to further constrain the planetary orstellar origin of the signals in addition to a classical periodogram analysis of an extensive number of activity indicators. We showed that with good statistical models, priors from auxiliary data, and elaborate simulations, planetary signals can be recovered and modeled that are on the order of the noise or even slightly weaker. Theproperties of all our detections are consistent with current formation scenarios regarding super-Earth planets.

Acknowledgments

We thank the anonymous referee, whose comments improved this work. This work was supported by the DFG Research Unit FOR2544 “Blue Planets around Red Stars”, project no. RE 2694/4-1. CARMENES is an instrument for the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, 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 Economy, 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. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia, Innovación y Universidades and the ERDF through projects PID2019-109522GB-C51/2/3/4 PGC2018-098153-B-C33 AYA2016-79425-C3-1/2/3-P, ESP2016-80435-C2-1-R, and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Canarias (SEV-2015-0548), Instituto de Astrofísica de Andalucía (SEV-2017-0709), and Centro de Astrobiología (MDM-2017-0737), the Generalitat de Catalunya/CERCA programme, and the NASA Grant NNX17AG24G. LCOGT observations were partially acquired via program number TAU2019A-002 of the Wise Observatory, Tel-Aviv University, Israel. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). The analysis of this work has made use of a wide variety of public available software packages that are not referenced in the manuscript: Exo-Striker (Trifonov 2019), astropy Astropy Collaboration (2018), scipy (Virtanen et al. 2020), numpy (Oliphant 2006), matplotlib (Hunter 2007), tqdm (da Costa-Luis 2019), pandas (The pandas development team 2020), and seaborn (Waskom et al. 2020).

Appendix A Priors for juliet

Appendix B Corner plots

Fig. B.1 Corner plot of the instrumental and GP parameters for GJ 251. Error bars denote the 68% posterior credibility intervals.

Fig. B.2 Corner plot of the planetary parameters for GJ 251. Error bars denote the 68% posterior credibility intervals.

Fig. B.3 Corner plot of the instrumental and GP parameters for HD 238090. Error bars denote the 68% posterior credibility intervals.

Fig. B.4 Corner plot of the planetary parameters for HD 238090. Error bars denote the 68% posterior credibility intervals.

Fig. B.5 Corner plot of the instrumental and GP parameters for Lalande 21185. Error bars denote the 68% posterior credibility intervals.

Fig. B.6 Corner plot of the planetary parameters for Lalande 21185. Error bars denote the 68% posterior credibility intervals.

Appendix C AliasFinder plot of HD 238090

Fig. C.1 Alias tests for HD 238090. The top plot (a) shows simulations motivated by a sampling frequency of . The bottom plot (b) shows simulations motivated by a sampling frequency of . Each row in these plots corresponds to one set of simulations for which the frequency of the injected signal is indicated by a vertical dashed blue line. The first row shows simulations with a period of 13.6838 d, and the second and third row show the simulations where the first-order aliases of 13.68 d, regarding the investigated sampling frequency were injected. Each column shows informative ranges of the periodograms, which are based on the assumed sampling frequency, and can be used for the comparison of data and simulations. From 1000 simulated data sets each, the median of the obtained periodograms (solid black line), the inter-quartile range and the ranges of 90% and 99% (gray shades) are shown. For comparison, the periodogram of the observed true data is plotted with a solid red line. The angular mean of the phase of some peaks and their standard deviation are shown in the clockdiagrams (black line and gray shades) and can be compared to the phase of these peaks in the observed periodogram (red line).

References

All Tables

All Figures

	Fig. 1 Posterior distribution in the αGP vs. Prot plane of the GP fit to the combined photometric data of T90, TJO, LCO, and SuperWASP for GJ 251. The color-coding shows the log-likelihood normalized to the highest achieved log-likelihood value within the posterior sample. Gray samples indicate solutions with Δ ln L lower than 10.
In the text

	Fig. 2 Joint GP model of the nightly binned photometric data of GJ 251. From top to bottom: T90 V, T90 R, TJO R, and SuperWASP.
In the text

Fig. 3 Generalized Lomb-Scargle periodograms of several activity indicators of GJ 251 from CARMENES spectroscopic data. The dashed black periodograms represent the GLS of the activity indicators, and the solid GLS periodogram represents the residuals after subtracting a 365 d signal. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The red dashed lines mark the rotation period and the first harmonic estimated from photometric data, while the dotted red lines show the 3σ uncertainties. The dashed black line marks a significant HIRES signal around 600 d, and the yearly period of 365 d. The dashed yellow line marks the period of the planetary signal published in this work.

In the text

Fig. 4 Generalized Lomb-Scargle periodograms of RV data for GJ 251 for CARMENES, HIRES, and a combination of both. The stellar rotational period derived by photometry is plotted as the dashed red line, and the 1σ and 3σ uncertainties are highlighted in red. We also indicate the harmonic of the rotational period and its uncertainty. The green line marks the suspected planetary signal. The blue lines mark the periods of the published planetary candidates by Butler et al. (2017) at 1.7 d, and the significant HIRES signal around 604 d.

In the text

Fig. 5 Stacked-Bayesian GLS periodogram of the planetary signal at 14.22 d and on the zero-planet residuals (left diagram with its ownnormalization). The three s-BGLS on the right were calculated using the one-planet Keplerian RV residuals and show signals that we attributed to stellar activity: the forest of signals between 60 and 75 d corresponding to roughly half of the rotation period (left), the RV signals around the photometrically derived rotation period at 119.5 d (middle), and the long period signal around 600 d. In all four plots, the dashed black line indicates the boundary between HIRES and CARMENES data, which were taken successively.

In the text

Fig. 6 Posterior distribution of the GP fit to the RV data in the αGP vs. Prot plane for GJ 251. The color-coding shows the log-likelihood normalized to the highest value in the posterior sample. Gray samples indicate solutions with a Δ ln L > 10 compared to the best solution. Top: GP fit to the RV data with a wide uniform prior for the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric GP results and additional constraints on the other hyperparameters. We overplot the derived rotational period of the photometric GP and its 3σ uncertainties with vertical lines. The horizontal line marks the cut in αGP used to constrain the GP fit shown in the lower plot.

In the text

Fig. 7 Gaussian process model for the RV data of GJ 251 without the planetary model, which is subtracted from the RV data. The constrained GP model is shown in red. The blue regions show 1σ, 2σ, and 3σ uncertainties. We show a zoom to some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP level of0.1%. In addition, we show the unconstrained GP model as the dashed black line in the upper plots and as the gray periodograms in the lower plots.

In the text

	Fig. 8 Left: radial velocity data with a combined model of one-planet and stellar activity using a Keplerian model and a GP. Right: plot phased to the orbital period of GJ 251 b without a GP component. The bottom plots show the residuals after the fit.
In the text

	Fig. 9 TESS PDCSAP light curve of GJ 251 folded to the 0.165 d signal. The black dots and uncertainty bars represent binned TESS data.
In the text

	Fig. 10 Posterior distribution of the GP fit to the photometric data in the αGP vs. Prot plane for HD 238090. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10 compared to the best solution.
In the text

	Fig. 11 Gaussian-process model for each nightly binned photometric data set of HD 238090. From top to bottom: NSVS, MEarth 2009 to 2010, MEarth 2010 to 2014, SNO V, SNO R, and TJO. For each instrument, we fit individual GP hyperparameters for the amplitudes and Γi, but we used global GP hyperparameters for the timescale of the amplitude modulation and the rotation period.
In the text

Fig. 12 Generalized Lomb-Scargle periodograms of several activity indicators based on spectroscopic data obtained by CARMENES for HD 238090. The dashed black periodograms represent the GLS on the activity indicators, and the solid GLS periodogram represents the residuals from which a 365 d sinusoidal signal was subtracted. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The solid red lines mark the rotation period estimated by photometric data, and the dashed red lines show the 3σ uncertainties. The dashed black line marks a period of 365 d, and the dashed yellow line marks the period of the planetary signal published in this work.

In the text

	Fig. 13 Generalized Lomb-Scargle periodogram of RV data of HD 238090 of the zero-planet fit, one-planet fit, and one-planet GP simultaneous fit. The stellar rotational period derived by photometry is plotted as the dashed red line, and 1σ and 3σ uncertainties are highlighted in red. The green line indicates the planet period at about 13.7 d.
In the text

	Fig. 14 Stacked-Bayesian GLS of the planetary signal (left) and a signal close to the estimated stellar rotation (right).
In the text

Fig. 15 Posterior distribution of the GP fit to the RV data in the αGP vs. Prot plane for HD 238090. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10. Top: GP fit to the RV data with a wide uniform prior to the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric GP results and an upper αGP constraint.

In the text

Fig. 16 Gaussian-process model for the RV data of HD 238090. The planetary signal is not included in the model and subtracted from the RV data. The constrained GP model is shown in red. The blue regions shows the 1σ, 2σ, and 3σ uncertainties. We show a zoom into some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP of 0.001. We also show the unconstrained GP model as the dashed black line in the upper plots and the gray periodograms in the lower plots.

In the text

	Fig. 17 Left: radial velocity data with a combined model for one planet and activity using a Keplerian model and a GP. Right: plot phased to the orbital period of HD 238090 b without a GP component.
In the text

Fig. 18 Generalized Lomb-Scargle periodograms of several activity indicators based on spectroscopic data obtained by CARMENES for Lalande 21185. The thinner dashed black periodograms represent the GLS on the activity indicators, and the solid GLS periodogram represents the residuals from which a 365 d sinusoidal signal was subtracted. For the residuals from which the 365 d signal was subtracted, we also overplot the s-BGLS periodogram, where the probability increases from blue to white to red. The dashed red lines mark the rotation period and the first harmonic estimated by photometric data, and the dotted red lines show the 3σ uncertainties. The dashed black lines mark a period of 365, 1400, and 2800 d, and the dashed yellow line marks the period of the planetary signal published in this work.

In the text

Fig. 19 Generalized Lomb-Scargle periodograms of Lalande 21185 for the different data sets and the combined CARMENES and SOPHIE data. Residual periodogram of the one-planet fit and the one-planet GP simultaneous fit for the combined CARMENES and SOPHIE data are also shown. Marked frequencies represent the claimed planet at 12.95 d by (Díaz et al. 2019; green), photometric stellar rotational period (red), the claimed planet candidate by (Butler et al. 2017; orange), and the long-term period (blue).

In the text

Fig. 20 Alias tests for Lalande 21185. The top plot (a) shows simulations motivated by a sampling frequency of . The bottom plot (b) shows simulations motivated by a sampling frequency of fs2 = 1.0027 d−1. Each row in these plots corresponds to one set of simulations for which the frequency of the injected signal is indicated by a vertical dashed blue line. The first row shows simulations with a period of 12.95 d, and the second and third row show the simulations in which the first-order aliases of 12.95 d, regarding the investigated sampling frequency, were injected. Each column shows informative ranges of the periodograms based on the assumed sampling frequency and can be used to compare data and simulations. From 1000 simulated data sets, we show the median of the obtained periodograms (solid black line), the interquartile range, and the ranges of 90% and 99% (gray shades). The periodogram of the observed data is plotted with a solid red line. The angular mean of the phase and the standard deviation is shown in the clock diagrams (black line and gray shades) and can be compared to the phase ofthe signals in the observed periodogram (red line).

In the text

Fig. 21 Stacked-Bayesian GLS periodogram of Lalande 21185. The color bar to the left color-codes each observation with the associated spectrograph (orange: CARMENES, and teal: SOPHIE). Left: s-BGLS of the zero-planet model at the period of theplanetary signal. Middle: s-BGLS on the residuals of the one-planet fit with a period of 12.95 d shown around a period of 55 d, which is the rotational period determined by Díaz et al. (2019). Right: s-BGLS on the residuals of the one-planet fit around the observed long-period signal.

In the text

Fig. 22 Posterior distribution in the GPα vs. Prot plane for Lalande 21185. The color-coding shows the log-likelihood normalized to the highest value within the posterior sample. Gray samples indicate solutions with a Δ ln L > 10. Top: GP fit to the RV data with a wide uniform prior to the rotational period. Bottom: GP fit to the RV data with an informative normal prior based on the photometric rotational period proposed by Díaz et al. (2019).

In the text

Fig. 23 Gaussian process model for the RV data of Lalande 21185. The planetary signal is not included in the model and subtracted from the RV data. The constrained GP model is shown in red. The blue regions shows the 1σ, 2σ, and 3σ uncertainties. We show a zoom into some CARMENES observations (top right). The GLS is evaluated on the GP model at each observed data point (top GLS) and daily (bottom GLS). The dashed line in the GLS periodograms indicates an FAP of 10−3. We also show the unconstrained GP model as the dashed black line in the upper plots and as the gray periodograms in the lower plots.

In the text

	Fig. 24 Left: radial velocity data with a combined model for one planet and activity using a Keplerian model and a GP. Right: plot phased to the orbital period of Lalande 21185 b without a GP component.
In the text

	Fig. 25 Confirmed exoplanets around M-dwarf host stars, detected with the RV method and with minimum masses below 10 M⊕ and orbital periods between 1 d and 400 d. The color-coding represents the spectral type. GJ 251 b, HD 238090 b, and Lalande 21185 b are marked in red. Data taken from http://exoplanet.eu
In the text

	Fig. B.1 Corner plot of the instrumental and GP parameters for GJ 251. Error bars denote the 68% posterior credibility intervals.
In the text

	Fig. B.2 Corner plot of the planetary parameters for GJ 251. Error bars denote the 68% posterior credibility intervals.
In the text

	Fig. B.3 Corner plot of the instrumental and GP parameters for HD 238090. Error bars denote the 68% posterior credibility intervals.
In the text

	Fig. B.4 Corner plot of the planetary parameters for HD 238090. Error bars denote the 68% posterior credibility intervals.
In the text

	Fig. B.5 Corner plot of the instrumental and GP parameters for Lalande 21185. Error bars denote the 68% posterior credibility intervals.
In the text

	Fig. B.6 Corner plot of the planetary parameters for Lalande 21185. Error bars denote the 68% posterior credibility intervals.
In the text

Fig. C.1 Alias tests for HD 238090. The top plot (a) shows simulations motivated by a sampling frequency of . The bottom plot (b) shows simulations motivated by a sampling frequency of . Each row in these plots corresponds to one set of simulations for which the frequency of the injected signal is indicated by a vertical dashed blue line. The first row shows simulations with a period of 13.6838 d, and the second and third row show the simulations where the first-order aliases of 13.68 d, regarding the investigated sampling frequency were injected. Each column shows informative ranges of the periodograms, which are based on the assumed sampling frequency, and can be used for the comparison of data and simulations. From 1000 simulated data sets each, the median of the obtained periodograms (solid black line), the inter-quartile range and the ranges of 90% and 99% (gray shades) are shown. For comparison, the periodogram of the observed true data is plotted with a solid red line. The angular mean of the phase of some peaks and their standard deviation are shown in the clockdiagrams (black line and gray shades) and can be compared to the phase of these peaks in the observed periodogram (red line).

In the text