© ESO 2020

1 Introduction

By July 2019, 665 multiplanetary systems were known, 148 of them discovered by precise Doppler spectroscopy. Only 11 of those Doppler spectroscopy-detected systems had stellar masses lower than 0.3 M_⊙. Such a low-mass star is YZ Ceti (GJ 54.1), which was reported to host three planets (Astudillo-Defru et al. 2017). It is the closest multiplanetary system to our Solar System published so far. Astudillo-Defru et al. (2017) announced that YZ Ceti is orbited by at least three Earth-mass planets at periods of 1.97 d, 3.06 d, and 4.66 d. The low amplitude of the signals (on the order of 1–2 m s⁻¹) together with the spectral window make the radial velocities of the system prone to strong aliasing. Although the system has been the subject of several studies (Astudillo-Defru et al. 2017; Robertson 2018; Pichierri et al. 2019) the true configuration of the planets is still highly disputed. For example, Robertson (2018) pointed out that the true configuration of the system could not be uniquely determined with the HARPS data available to Astudillo-Defru et al. (2017). In particular, the signal of planet c (P = 3.06 d) had a strong alias at 0.75 days.

Astudillo-Defru et al. (2017) also mentioned a fourth tentative signal slightly above a one-day periodicity at 1.04 d. On the other hand, Tuomi et al. (2019), who used only 114 radial velocity measurements by HARPS and 21 from HIRES, supported the existence of only two planet candidates for YZ Ceti at periods of 3.06 d and 4.66 d.

Determining the true configuration of a planetary system and constraining its parameters is of the utmost importance for any attempt to perform a dynamical characterization or to understand its formation. Recent planet formation results suggest that ultra-compact planetary systems around stars similar to YZ Ceti should be common, where the planets are typically locked in long resonant chains, exhibiting both two-body and three-body resonances (Coleman et al. 2019). We took additional radial velocity (RV) data for YZ Ceti with CARMENES and HARPS to address the open questions regarding the exact planetary configuration, the possibility of additional companions, the modeling, the influence of stellar activity, and the dynamical properties of this multiplanetary system.

This work is organized as follows. The data and instruments used in this study are described in Sect. 2. We discuss the basic stellar properties of YZ Ceti and the analysis of the photometry and activityindicators in Sect. 3. In Sect. 4 we examine the RV data for additional planet candidates, resolve the alias issues raised in the literature, and present strong evidence for the correct planetary configuration of the system by resolving the alias issues raised in the literature. In Sect. 5 we describe the modeling with a Gaussian process (GP) and the influence of the stellar activity on the eccentricity of the innermost planet, while in Sect. 6 we constrain our posterior parameters even further by adopting the criterion of long-term stability for the multiplanetary system. The results of this study are summarized and discussed in Sect. 7.

2 Data

2.1 CARMENES

We observed YZ Ceti as one of the 324 stars within our CARMENES¹ Guaranteed Time Observation program to search for exoplanets around M dwarfs (Reiners et al. 2018). CARMENES is a precise échelle spectrograph mounted at the 3.5 m telescope at the Calar Alto Observatoryin Spain. It consists of two channels: the visual (VIS) covers the spectral range 0.52–0.96 μm with spectral resolution of R = 93 400, and the near-infrared (NIR) the 0.96–1.71 μm range with spectral resolution of R = 81 800 (Quirrenbach et al. 2014). We obtained 111 high-resolution spectra in the VIS and 97 spectra in the NIR between January 2016 and January 2019. Three spectra of the VIS and NIR arm were without simultaneous Fabry-Pérot drift measurements and were therefore excluded from our RV analysis. The data were reduced using CARACAL (Caballero et al. 2016b), and we obtained the radial velocities using SERVAL (Zechmeister et al. 2018). SERVAL determines RVs by coadding all available spectra of the target with signal-to-noise ratio (S/N) higher than 10 to create a high S/N template of the star and by deriving the relative RVs with respect to this template using least-squares fitting. The radial velocities were corrected for barycentric motion, secular perspective acceleration, instrumental drift, and nightly zero points (see Trifonov et al. 2018; Tal-Or et al. 2019, for details). For the VIS we achieved a median internal uncertainty, including the correction for nightly zero points, of 1.55 m s⁻¹ and a root mean square (rms) of 3.00 m s⁻¹ around the mean value. For the NIR we achieved a median internal uncertainty, including the correction for nightly zero points, of 5.71 m s⁻¹ and an rms of 7.17 m s⁻¹ around the mean value. Due to the small amplitudes of the planetary signals of less than 2 m s⁻¹ in the YZ Ceti system, we only used the VIS data for our analysis.

The radial velocity time series and their uncertainties for all data sets used within this work are listed in Table B.1.

2.2 HARPS

The High Accuracy Radial velocity Planet Searcher (HARPS) (Mayor et al. 2003) is a precise échelle spectrograph with a spectral resolution of R = 110 000 installed at the ESO 3.6 m telescope at La Silla Observatory, Chile. HARPS covers the optical wavelength regime, and was the first spectrograph that reached a sub-m s⁻¹ precision. We retrieved 334 high-resolution spectra from the ESO public archive, of which 59 were collected by the Red Dots program (Dreizler et al. 2020). As with the CARMENES data we used SERVAL to obtain the RV from the corresponding spectra. Only 326 spectra were used for the coadding of the template, as eight spectra had a S/N of less than 10. We then calculated the corresponding RVs for all 334 spectra. From these 334 we removed two extra measurements: one at BJD 2 456 923.73068 as it was an obvious outlier with very low S/N of 3.6 and one at BJD 2 458 377.92388 due to an RV uncertainty larger than 83 ms⁻¹. This resulted in atotal of 332 RV measurements by HARPS. We divided the HARPS RV data due to a fiber upgrade on May 28, 2015 (Lo Curto et al. 2015), into pre- and post-fiber data and fitted an offset for it. For the pre-fiber upgrade data sets we achieved a median internal uncertainty of 1.92 m s⁻¹ and an rms of 2.88 m s⁻¹ around the mean value. For the post-fiber upgrade data sets we achieved a median internal uncertainty of 1.86 m s⁻¹ and an rms of 3.72 m s⁻¹ around the mean value.

2.3 Photometry

Details of observations from five photometric facilties are given below.

ASAS

The All-Sky Automated Survey (Pojmański 1997) has been monitoring the entire southern and part of the northern sky since 2000. We retrieved 461 ASAS-3 measurements of YZ Ceti taken between November 2000 and November 2009.

ASAS-SN

The All-Sky Automated Survey for Supernovae (ASAS-SN) (Shappee et al. 2014; Kochanek et al. 2017) uses 14 cm aperture Nikon telephoto lenses, each equipped with a 2048 × 2048 ProLine CCD camera with a field of view (FOV) of 4.5 deg² at different observatories worldwide. We extracted about six years of photometric observations (2013–2019) in the V band from the ASAS-SN archive² for YZ Ceti.

MEarth

The MEarth project (Berta et al. 2012) is an all-sky transit survey. Conducted since 2008 it uses 16 robotic 40 cm telescopes, 8 located in the northern hemisphere at the Fred Lawrence Whipple Observatory in Arizona, USA, and 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 2048 × 2048 CCD that provides a FOV of 26 arcmin². MEarth generally uses an RG715³ long-pass filter, except for the 2010-2011 season when an I_715−895 interference filter was chosen. In the case of YZ Ceti, we used archival data from MEarth telescopes T11 and T12 released in the seventh data release (DR7⁴). The set from T11 consists of 40 epochs with a time span of 102 days between May and August 2017, while the set from T12 consists of 25 epochs with a time span of 68 days between June and August 2017.

ASH2

The Astrograph for Southern Hemisphere II (ASH2) telescope is a 40 cm robotic telescope located at the San Pedro de Atacama Celestial Explorations Observatory (SPACEOBS) in Chile. The telescope is equipped with a STL11000 2.7k ×4k CCD camara and has a FOV of 54 × 82 arcmin². We carried out the observations in the V and R bands and in two runs. Run 1 consisted of 14 observing nights during the period September to December 2016, with a time span of 65 days and about 650 data points collected in each filter. Run 2 consisted of 28 observing nights between July and October 2018, with a time span of 91 days and about 500 data points collected in each filter.

SNO

The Sierra Nevada Observatory (SNO) in Spain operates four telescopes. The T90 telescope at Sierra Nevada Observatory is a 90 cm Ritchey-Chrétien telescope equipped with a CCD camera VersArray 2k × 2k, FOV 13.2 × 13.2 arcmin (Rodríguez et al. 2010). We carried out the observations in both Johnson V and R filters on 38 nights during the period August to December 2017, with a time span of 125 days and about 1850 data points collected in each filter.

3 Properties of YZ Ceti

3.1 Basic physical parameters

YZ Ceti (GJ54.1) is an M4.5 V star at a distance of approximately 3.7 pc (Gaia Collaboration 2018), making it the 21st nearest star to the Sun⁵. We provide an overview of the basic stellar parameters in Table 1. For our analysis we adopted the stellar parameters of Schweitzer et al. (2019). The effective temperature T_eff, surface gravitiy log g, and metallicity [Fe/H] were determined by fitting PHOENIX synthetic spectra (Husser et al. 2013) to CARMENES spectra using the method of Passegger et al. (2018). The luminosity was derived by integrating broadband photometry and adopting the parallax measurement from Gaia DR2. The stellar radius R of 0.157 ± 0.005 R_⊙ was determined by applying the Stefan-Boltzmann law. With a linear mass-radius relation we then obtained 0.142 ± 0.010M_⊙ for the stellar mass (see Schweitzer et al. 2019, for details). We computed the Galactocentric space velocities UV W as in Cortés-Contreras (2016) with the latest equatorial coordinates, proper motion, and parallax from Gaia DR2 and the latest absolute RV of Lafarga et al. (2020). With these UV W values, YZ Ceti kinematically belongs to the Galactic thin disk and has never been assigned to a stellar kinematic group (Cortés-Contreras 2016). Apart from flaring activity, which is very frequent in intermediate M dwarfs of moderate ages, the star does not display any feature of youth. Engle & Guinan (2017) estimated an age of 3.8 ± 0.5 Gyr from the stellar rotation and X-ray emission, which is consistent with the star kinematics.

3.2 Photometric analysis

Astudillo-Defru et al. (2017) estimated the rotation period of YZ Ceti to be 83 d by analyzing ASAS photometry and the FWHM of the cross-correlation function (FWHM_CCF) of the radial velocity data obtained by HARPS. However, shortly after publication Jayasinghe et al. (2017) determined the photometric rotation period at P_rot = 68.3 d using an ASAS-SN lightcurve of 854 photometric measurements in the V band obtained between 2013 and 2017. This result was confirmed by Engle & Guinan (2017), who estimated a rotation period of P_rot = 67 ± 1.8 d based on V -band photometry taken between 2010 and 2016 with the 1.3 m Robotically Controlled Telescope. Both of these estimates are in good agreement with the value of P_rot = 69.1 d determined by Suárez Mascare no et al. (2016) and P_rot = 69.2 ± 0.4 d determined by Díez Alonso et al. (2019). Furthermore, Fig. 1 in Astudillo-Defru et al. (2017) shows that the second highest peak of the periodogram of the ASAS data, as well as the highest peak of the FWHM_CCF between JD 2 457 100 to JD 24 577 000, was around 68 days.

We combined the public archive data from ASAS, ASAS-SN, and MEarth telescopes 11 and 12 (Sect. 3.2) with our own observationswith ASH2 and SNO T90 to carry out the most extensive combined photometric analysis of the rotation period of YZ Ceti. For each instrument and photometric filter we created generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) on the nightly binned data. We show the obtained GLS peaks for each instrument in Table 2. All instruments except ASAS directly show a highly significant peak at around 68 d. However, for some instruments the highest peak was either a yearly alias or close to half the rotational period of this ~68 d periodicity. The formal uncertainties on the frequency, and therefore for the period for each peak, are estimated by the GLS routine from the local χ² curvature. This estimate does not account for incorrect choices of alias peaks, hence real uncertainties can be larger. We have also done a combined analysis of all instruments in the R and V bands by fitting for an offset and jitter term for each instrument within the two different bands. We also do a combined analysis for the SNO R-band data and the MEarth T 11 as they have no peak at 68 d; instead, they have a single broad peak close to the yearly alias, which however includes the 68 d period. Combining both data sets yields 41.63 ± 0.85 d as the highest and 69.31 ± 2.52 d as the second highest peak, both with a false-alarm probability (FAP) below 0.001.

In Fig. 1 we show the nightly binned photometric time series of each instrument in the V band and R band as well as the phase plots corresponding to each time series. The bottom row of panels displays the GLS periodograms of each single instrument in the corresponding photometric band (left V, middle R and I) and a combination of all instruments in the V band and combined R and I bands.

In some of our data sets we recognized periodicities near 80 d and 57 d. The former period is close to the rotation period adopted by Astudillo-Defru et al. (2017). Fitting a sinusoid for either the 68-day signal or the 80-day signal removes the other, a strong sign of aliasing. In all data sets the 80-day signal can be reproduced by strong aliases due to annual sampling effects in the window function together with the 68-day periodicity. There are no strong signs of the 80-day periodicity in the combined R-band data nor in many individual R- and V -band data sets. These subsets did not show large annual peaks in the window function, so we are confident that the true rotation period is around 68 days.

From the R + I band combined GLS analysis we determined a rotation period of 68.5 ± 1.00 d with an amplitude of 8.6 ± 0.7 mmag, and we independently determined from the V band a rotation period of 68.40 ± 0.05 d with an amplitude of 14.2 ± 0.5 mmag.

3.3 Spectroscopic activity indicators

YZ Ceti is an active M dwarf identified as a flare star (Kunkel 1970; Shakhovskaya 1995; Montes et al. 2001). Reiners et al. (2018) estimated an upper limit of v sin i < 2 km s⁻¹, which corresponds to a slow rotational velocity. The equatorial rotation speed estimated from the radius and rotation period of the star is 2πR sin i∕P ≈ 0.12sin(i) km s⁻¹, well below the directly estimated upper limit of 2 km s⁻¹. In addition to the photometric observations, we analyzed several spectral activity indicators from the CARMENES and HARPS spectra. We searched for periodicities of the chromatic index (CRX), Hα, and differential line width (dLW) within all data sets (see Zechmeister et al. 2018, for CRX and dLW), and the Ca II IRT lines within the CARMENES data. These indicators are directly provided by SERVAL. From the CARMENES data we also determined the full width half at maximum (FWHM), contrast, and bisector span (BVS) of the cross-correlation function (CCF) (see Reiners et al. 2018). The GLS periodograms for each of these indicators are shown in Fig. 2. We identified a forest of significant peaks that include both the 80 d and 68 d periods visible from photometry within the CARMENES CRX, CARMENES dLW, CARMENES FWHM, CARMENES contrast. Within the HARPS activity indicators we identify a significant peak in the HARPS-POST dLW at 69.68 ± 0.23 d, where the error represents the 1σ uncertainty. These indicators (CRX, dLW, and FWHM) are sensitive to the photosphere of the star and are in agreement with our derived photometric rotation period of the star. However, we did not find a significant correlation between the CARMENES RVs and CARMENES CRX at this period. We did not see any significant signals for the remaining activity indicators. The forest of significant signals could be explained by adopting a period of around 68 days and calculating possible yearly aliases due to the window function of the radial velocity observations. Overall, we found a good agreement between the spectral activity indicators and the derived photometric rotation period. In particular, we did not find any signs of activity close to the periods of the planetary signals, or their aliases (see below).

4 YZ Ceti planetary system

4.1 Search for planetary signals

Compared to the discovery paper with 211 data points, we obtained 121 additional data points from HARPS and 108 from CARMENES, resulting in a total of 440 radial velocities for YZ Ceti, which more than doubles the number of available RV measurements in previous studies of this system; for example, both Astudillo-Defru et al. (2017) and Robertson (2018) used only 211 data points taken by HARPS before October 2016. For a major fraction of the CARMENES observations, we took two observationsper night in order to break possible degeneracies due to daily aliasing. In Fig. 3 we show the RV data and the final stable fit (see Sect. 6). For the fitting of planetary signals, we used the tool juliet (Espinoza et al. 2019), which allows fitting of photometric and/or RV data by searching for the global maximum of the Bayesian evidence within the provided prior volume of the fitting parameters. It does so by using nested sampling algorithms, for example MultiNest (Feroz et al. 2009), PyMultiNest (Buchner et al. 2014), and dynesty (Speagle 2020). For the modeling, juliet usesmany different publicly available packages, for example batman (Kreidberg 2015) for transits, radvel (Fulton et al. 2018) for radial velocities, and george (Ambikasaran et al. 2015) and celerite (Foreman-Mackey et al. 2017)for GP. The priors of the planetary parameters for the fits are shown in Table C.1.

We did a periodogram analysis of the RV data and fitted for the strongest signal until no significant peak with a FAP of less than 0.001 was observed in the residuals periodogram. The FAPs were determined by bootstrapping using 10 000 realizations. In order to test whether fitting n + 1 planets was significantly better than fitting for n planets, we compare the Bayesian log-evidence of the corresponding models. The resulting log-evidence for each model are displayed in Table 3. As suggested by Trotta (2008), we regard the difference between two models as strongly significant if their log-evidence differs by $\Delta \ln \mathcal{Z}>5$. The residual periodograms for these runs are shown in Fig. 4; the strongest signal in the periodogram was at 3.06 d. After fitting this signal, the next highest peak was at 4.66 d, and thereafter at 2.02 d. After fitting for three planets, the resulting periodogram did not show any remaining significant signal (no signal with FAP < 0.01; see also Fig. 4). Examining the orbital parameters of the three planets, we found that planet b at 2.02 d has an unusually high eccentricity $e_b={0.41}_{-0.17}^{+0.14}$ [0.12, 0.67], where the errors are 1σ uncertainties and the values inside the bracket represent the 95% density interval. The eccentricity is not significantly different when we choose the alias at 1.97 d, which was the favored period by Astudillo-Defru et al. (2017) for YZ Ceti b. The high eccentricity of most of the posterior samples led to instability for the majority of the samples on very short timescales as shown by integrations using an N-body integrator (see Sect. 5 and Table 7). Therefore, we compared the Bayesian log-evidence of different configurations where we fixed the eccentricity for all or certain combinations of planets to zero (shown in Table 3). This procedure reduced the number of parameters per planet from five to three (as ω is not defined for e = 0). We find that a model where we fix the eccentricity for planet c to zero but keep e_b and e_d open performs moderately better than a fit that fits eccentric orbits for all planets ($\Delta \ln \mathcal{Z}=3>2.5$) or fixes eccentricity to zero for all planets ($\Delta \ln \mathcal{Z}=2.7>2.5$). As a result, there is moderate evidence to fit eccentric orbits for planets b and d (Trotta 2008). The residual periodogram of a three-planet circular fit showed several remaining peaks. The strongest of them was at 67.69 d, and then the peaks at 69.22 d and 68.28 d, all with FAP < 0.01, and at 0.98 d and 1.02 d with a FAP < 0.1. All these peaks could be directly attributed to the stellar rotational period, which is known from photometry to be 68.4 d. The peaks at 67.69 d and 69.22 d are yearly aliases of the 68.28 d period, and 1.02 d and 0.98 d are its daily aliases. We also identified a peak at 28.37 d with FAP < 0.1. However, by fitting a simple sinusoid to the activity signal around 68 d, the 28.37 d singal is removed, a strong indication that this signal is connected to activity.

We tested the coherence of the 28.37 d signal over time with the stacked-Bayesian GLS periodogram (s-BGLS) method (Mortier et al. 2015; Mortier & Collier Cameron 2017). These BGLS periodograms allow comparison of the probabilities of the signals with each other, while the stacking allows assessment of the coherence of the signal with increasing number of observations. As in Mortier & Collier Cameron (2017) we normalized all s-BGLS periodograms to their respective minimum values.

We found that the 28.37-day signal (f ≈ 0.035 d⁻¹) was not very stable over the observed time interval, which can be seen in Fig. 5. In particular, this period was not prominent within CARMENES spectra. From these BGLS periodogograms we also deduce that the activity related to the 68-day signal increased over time.

From our GLS analysis we also identified a signal at 7.05 d with FAP just slightly above a FAP of 0.1 in the residuals of the three-planet circular fit. After fitting a three-planet model simultaneously with a GP to model the activity (see Sect. 5), we find that this signal increases slightly; it is the highest remaining signal in the residuals, but still below our detection threshold. The signal has an amplitude of roughly 0.6 m s⁻¹, and is close to an optimal 3:2 commensurability with regard to the period of YZ Ceti d, hinting at the possibility that there might be a fourth planet in the system. However, the signal requires more data to either confirm or refute its presence.

With the data analyzed to date we thus cannot find statistically significant evidence for an additional fourth planet in the YZ Ceti system. In particular, the tentative signal mentioned by Astudillo-Defru et al. (2017) at 1.04 d (f ≈ 0.962 d⁻¹) has decreased in significance (see also Fig. 5, left). In contrast to these signals, our s-BGLS analysis shows that the additional observations increase the probability of the signals at the frequencies of the three planets, further strengthening the possibility of a planetary origin of these signals.

Fig. 1 Top three rows: nightly binned photometric time series (top: all data in the V band, middle: zoom to ASH2 and SNO data in the V band, bottom: all data in the R band) and aphase plot to the determined rotation period. Bottom row: periodograms of the different instruments in the V band (left), R band (middle), and the combination of all instruments in each band (right). The periodograms for the analysis on individual instruments are color-coded (blue: ASH2, brown: SNO, red: ASAS, green: ASAS-SN, dark green: MEarth T11, purple:MEarth T12), while the combined periodograms are plotted in black. The solid line represents the combined V -band periodogramand the dashed line the combined R-band periodogram. For the combined periodograms we show the FAP level of 0.001 (green solid line for V band and red dashed line for R band). The vertical black line in each periodogram represents the determined rotation period of 68.4 d and 68.5 d, respectively.

Fig. 2 GLS periodograms of several activity indicators in the CARMENES and HARPS data. The periodograms are separated and show different frequency regimes to better display the significant peaks within the low-frequency regime. The red dashed line shows the photometric rotation period and the black solid lines highlight the planetary frequencies.

4.2 Configuration of the system

Each of the three planetary signals had at least one strong alias, making it difficult to pin down the correct period for the planets. The three strong alias pairs are P_b = 1.97∕2.02 d, P_c = 0.75∕3.06 d, and P_d = 1.27∕4.66 d. Astudillo-Defru et al. (2017) published the configuration P_b = 1.97 d, P_c = 3.06 d, and P_d = 4.66 d while Robertson (2018) favored the orbital configuration P_b = 1.97 d, P_c = 0.75 d and P_d = 4.66 d. Here we show that the most likely configuration of YZ Ceti is P_b = 2.02 d, P_c = 3.06 d, and P_d = 4.66 d and that we can robustly determine the configuration of the system with our data. We use two distinct methods. First, we compared the maximum log-likelihood, similar to Robertson (2018), of different realizations of periods for the three-planet models by sampling each possible configuration within juliet. When we fitted for the different aliases we only changed the prior of the period and ensured that the posterior is well sampled and not truncated within this volume; for example, to fit planet c at 0.75 d instead of 3.06 d we simply adopted $\mathcal{U}(0.7,0.8)$ instead of $\mathcal{U}(3.0,3.1)$. The log-likelihood is more robust to changes in the prior than the log-evidence when comparing the same model with equal number of free parameters but different realizations. The achieved maximum log-likelihoods for each sample of this analysis are summarized in Table 4. The data significantly favors a model with planets at 2.02 d, 3.06 d, and 4.66 d ($\Delta \ln \mathcal{L}=6.5>5$, with respect to the best model). In particular, the models with the proposed alias by Robertson (2018) adopting P_c = 0.75 d perform worse.

As outlined by Dawson & Fabrycky (2010), the best fitting model does not necessarily represent the true configuration of the system. Therefore, model comparison can only be a strong indication of the way to disentangle aliases. Hence, we also applied a slightly modified version of the method described by Dawson & Fabrycky (2010). The basic idea is to simulate periodograms of the candidate frequencies and their possible aliases by investigating the periodogram of the spectral window function (Roberts et al. 1987), which we show for our data in Fig. 6. Using the same sampling as for the original data, a signal with the same properties (phase, amplitude period) is injected into a simulated time series. A periodogram analysis is used to compare the signal properties at the proposed true frequency and at each alias frequency of the simulated and the data periodograms. If the periodogram of one of the simulated frequencies matches the observed data significantly better than the others, this frequency can be considered to be the true one. The original Dawson & Fabrycky (2010) method did not include any noise for the simulated periodograms. Dawson & Fabrycky (2010) stated that only if the noiseless simulated periodograms matches the data periodogram well can it be regarded as a good match. Such a noiseless periodogram, which only includes the true frequency, tends to look very clean with sharp peaks. This makes it in some cases difficult to compare with the noisier data periodogram. In an attempt to break the period degeneracies for YZ Ceti, Robertson (2018) used the method by Dawson & Fabrycky (2010) and included noise by adopting the uncertainties of the radial velocities at each time together with a white noise model of the star based on its derived jitter value. Nevertheless, with the data available at that time, Robertson (2018) was unable to constrain the true frequencies of YZ Ceti b and c with this test. Care should always be taken if noise is included in the simulations. Dawson & Fabrycky (2010) already mentioned that noise can interfere with the candidate periods resulting in the possibility that the power of the alias frequency is higher than the power of the true frequency inthe data periodogram. Creating only one realization, a simulated periodogram with noise can therefore lead to incorrect conclusions. Most importantly, however, the inclusion of noise can have a significant effect on the derived phases, as small errors of the determined phase value can accumulate to significant phase differences for peaks that are far away from the injected signal in frequency space. This can be easily seen and tested by injecting two signals with slightly different phases into a simulated time series and examining the phase of their aliases.

Therefore, if noise is included it is necessary to also account for the uncertainties of the determined phase values. We recommend the addition of noise via this method, but suggest the following approach: coupling the method with a Monte Carlo technique and creating 1000 different versions of the simulated data sets. For each time series we used a white noise model, as in Robertson (2018), so that we draw for each realization i from a Gaussian distribution with ${\sigma }_i^2={\sigma }_{\text{RV},i}^2+{\sigma }_{\text{jitter}}^2$. To compare the simulated data with the observations, we created a master periodogram, which is the median of the periodograms from all simulations. This was repeated for the expected true frequency and its first-order aliases. The aliases were thereby calculated using the equation $$f_{\text{alias}}=f_{\text{p}}\pm m\cdot f_{\text{s}},$$(1)

where f_p is the planetary orbital frequency, f_s the sampling frequency (in our case the largest peak of the window function periodogram), and m the order of the alias. We show the results from this analysis in Fig. 7.

The first row in each plot corresponds to the simulation for the strongest peak in the observed periodogram (the expected true frequency f_p), while the second and third row correspond to its first-order daily aliases of f_alias = f_p − f_s and f_alias = f_p + f_s respectively (f_s = 1.0027 d⁻¹). The signal that is most likely underlying the observed periodogram is the one whose simulated periodograms fits best all three subsets compared to the simulations in the other rows. Such a plot is recommended for any system where possible strong aliases exist. We provide the script used for the alias-testing on github⁶ (Stock & Kemmer 2020).

The key frequencies and periods for the following analysis are summarized in Table 5. We started with the strongest signal in our data at 3.06 d (f ≈ 0.327 d⁻¹) which was the signal called into question by Robertson (2018). The top panel of Fig. 7a shows a significantly better agreement between the simulated periodogram and the observed data when the 3.06-day signal is injected into our simulation compared to when we use its aliases. Injecting the alias at 1.48 d (f ≈ 0.676 d⁻¹, see Fig. 7a, middle panel) did reproduce the phase values of all peaks well, but not the peak height at 0.75 d (f ≈ 1.332 d⁻¹). Injecting the 0.75-day periodicity (Fig. 7a, lower panel) resulted in large phase differences for the other two aliases, and the peak height of the signal at 1.48 d was not well reproduced. From this result and the previousresults based on the Bayesian evidence we concluded that the true period of YZ Ceti c is indeed 3.06 d.

We subtracted this signal from the periodogram using a simple sinusoid and tested the candidate periods of YZ Ceti d, the second strongest signal in the data, for its possible aliases. In the top panel of Fig. 7b, the simulated periodogram with a period of 4.66 d (f ≈ 0.215 d⁻¹) fits the peak heights and phases of the data periodogram well. Injecting the signal of 1.27 d (f ≈ 0.788 d⁻¹; see Fig. 7b, middle panel) did not reproduce the peak at 0.82 days (f ≈ 1.218 d⁻¹) and resultedin larger differences in the phase values of the other two peaks. Instead, the 0.82-day signal (Fig. 7b, lower panel) did not reproduce the 1.27-day peak. We concluded that YZ Ceti d orbits at 4.66 d as stated by Astudillo-Defru et al. (2017) and Robertson (2018). We subtracted this signal and analyzed the aliases for YZ Ceti b.

Regarding the simulated periodograms for YZ Ceti b in the top panel of Fig. 7c, we found that the simulated periodogram with the period of 1.97 d (f ≈ 0.508 d⁻¹) published by Astudillo-Defru et al. (2017) performs worse than the alternative period of 2.02 d (f ≈ 0.495 d⁻¹; see Fig. 7c, middle panel). The 1.97-day signal did not reproduce the peak at 0.67 d (f ≈ 1.497 d⁻¹) at all, as the peak heights of the simulated periodograms deviated significantly from the data periodogram, and the data periodogram is not in the range of 99% of the simulated periodograms. The same is true for the peak at 1.97 d when simulating the 0.67-daysignal (Fig. 7c, lower panel). We therefore adopted a period of 2.02 d for YZ Ceti b, in line with our results by juliet and in contrast to the published period by Astudillo-Defru et al. (2017) and Robertson (2018). We also tested subtracting the alternative alias solutions from the periodograms before doing the analysis for YZ Ceti b and d. This had no significant influence on the results presented in Fig. 7.

With more data, some RV measurements taken twice per night and with observations taken from multiple observatories (Calar Alto, Spain and La Silla, Chile), which have a ~ 4 h difference inlongitude, we have improved sampling and therefore were able to solve the alias problem for this system. Overall, we found from our analysis on aliases and the model comparison within the framework of Bayesian evidence that the most probable configuration of the YZ Ceti multiplanetary system, as derived from the current data, is a three-planet system with planets at periods of 2.02 d, 3.06 d, and 4.66 d. Coincidentally, attributing YZ Ceti b to a period of 2.02 d instead of 1.97 d brings the system configuration closer to a 3:2 period commensurability for both pairs of neighboring planets.

Fig. 3 Radial velocity data and final stable fit including the Gaussian process model for the activity signal (see Sects. 5 and 6).

Fig. 4 Generalized Lomb-Scargle periodograms of the data and different fit models from Table C.1, and the final GLS periodogram after fitting for the activity with a GP (see Sect. 5).

Fig. 5 S-BGLS periodograms after subtracting the three planetary signals. Left: around the frequency of the tentative signal mentioned by Astudillo-Defru et al. (2017), middle: 29.36-day signal visible as remnant power in a circular three-planetmodel, right: 68-day signal attributed to the stellar rotation.

Fig. 6 Periodogram of the spectral window function for the data used in our analysis.

Fig. 7 Alias tests for the periods of 3.06 d (a), 4.66 d (b), and 2.02 d (c). In each plot,each row corresponds to one set of simulations. The frequency of the injected signal is indicated by a vertical blue dashedline. From 1000 simulated data sets each, the median of the obtained periodograms (black solid line), the interquartile range, and the ranges of 90 and 99% (shades of gray) are shown. For comparison, the periodogram of the observed data is plotted as a red solid line. Additionally, the angular mean of the phase of each peak and its standard deviation of the simulated periodograms are shown as clock diagrams (black line and grays) and can be compared to the phase of thepeak in the observed periodogram (red line).

5 Simultaneous fitting of stellar activity and Keplerians

As in Astudillo-Defru et al. (2017), we simultaneously modeled the stellar activity by using GP regression models as we found that correlated noise seems to influence the derived planetary parameters. This is in contrast to Robertson (2018), who did not use a GP to model the YZ Ceti system. Compared to sinusoidal signals, a GP has the advantage that is more flexible, and can therefore capture more features resulting from the stellar activity. However, GPs may also potentially lead to overfitting the data and so to absorbing noise or planetary signals into the presumed stellar activity. Weused juliet to model the activity signal simultaneously with the Keplerian models. We used an exp-sin-squared kernel multiplied with a squared-exponential kernel (Ambikasaran et al. 2015). This kernel has the form $$k_{i,j}={\sigma }_{{\text{GP}}_i}^2\exp (-{\alpha }_i{\tau }^2-{\Gamma }_i{\sin }^2(\pi \tau /P_{\text{rot},i})),$$(2)

where ${\sigma }_{{\text{GP}}_i}$ is the amplitude of the GP component given in m s⁻¹, α_i is the inverse length scale of the GP exponential component given in d⁻², Γ_i is the amplitude of GP sine-squared component given in m s⁻¹, P_rot,i the period of the GP quasi-periodic component given in d, and τ is the time-lag. The α value is a measure of the strength of the exponential decay of the quasi-periodic kernel. A lower α describes a more stable periodic signal and stronger correlation between data points. The quasi-periodic kernel is a kernel widely used in the literature for the modeling of stellar activity with a GP. We also tested the celerite approximation to a quasi-periodic kernel and the celerite Simple Harmonic Oscillator (SHO) kernel (Foreman-Mackey et al. 2017). Both of these had significantly worse evidence than the quasi-periodic kernel and were for most of the runs also not able to reproduce the rotational period of YZ Ceti.

We performed runs with different priors for the GP hyperparameters in order to investigate the influence of the GP modeling on the planetary parameters. YZ Ceti is a rather compact system, so we also used the dynamical stability of the derived orbits to test whether our posterior distributions for the planetary parameters are realistic (see Sect. 6). The following runs were all done by using the dynesty (Speagle 2020) dynamic nested sampling with 1500 live points and the sampling option slice, which is needed for our high-dimensional parameter space.

Our very first GP model used uninformative priors of the GP hyperparameters that spanned a very wide parameter range, especially for the rotation period, for which we used $\mathcal{U}(30,100)$. From the posterior samples of this run, we found that only for rotation periods around 68 d the GP model allowed very low α values, consistent with a rather stable periodic signal. This period range is consistent with the photometric rotation period. However, we also observed a plateau of a large number of possible solutions that range over the complete prior volume but have rather high values of α between 10⁰ d⁻² and 10⁻⁴ d⁻². Therefore, only a few samples modeled the stellar rotational signal.

In the four GP runs, which we describe in detail below, we constrained the prior on the rotation period to sample more dense around the period range of the photometric rotation period. The reason was that we wanted to use the GP primarily as a model for the stellar rotational signal and not for any other residual noise (e.g., instrumental). Our analysis showed thatthe correlated signal originating from the stellar rotational signal affects the planetary parameters the most, especially for YZ Ceti b at 2.02 d. We therefore note that the following approach might be unique for YZ Cetiand systems that suffer from a similiar problem. The GP priors used for these runs are listed in Table 6, where we also show the posteriors of the hyperparameters as well as the evidence of the run, the maximum achieved log-likelihood, and the eccentricity of planet b as this parameter is rather sensitive to the modeling of the stellar activity. The Bayesian log-evidence of all these runs was significantly better than a simple three-planet Keplerian fit to the data. In Fig. 8 we show the scatter plots of the sampled α values of the quasi-periodic kernel over the sampled rotational periods. Since the influence of the activity on the RV data is wavelength dependent, and HARPS and CARMENES operate across different wavelength regimes, we also tested for each run whether it is justified to use distinct GP amplitudes (σ_GP and Γ) for the two spectrographs separately. For all our runs with distinct GP amplitudes we achieved less log-evidence while increasing the number of parameters. The derived planetary parameters and remaining non-amplitude GP parameters were not significantly different from runs where we did not use distinct GP amplitudes for the instruments. Therefore, we stayed with the simplest GP model which has fewer hyperparameters and a higher log-evidence, but global amplitude hyperparameters.

For run a we set up a narrow uninformative uniform prior around the region of the suspected stellar rotational period. From the results of the posterior samples we found that the GP still did not predominantly model the rotation period of the star in most of the samples. The posterior of the GP rotational period was mostly flat and not well constrained, as [65.24, 70.00] which populates almost the complete provided prior range of this parameter. In Fig. 8a we see the same plateau as before, in the range of α values from ~ 0.1 d⁻² to 1 d⁻², corresponding to decay timescales of several days ($\tau ~{\alpha }_{\text{GP}}^{-1/2}$). The plateau spanned the range between 10⁰ d⁻² and 10⁻⁴ d⁻² in α for this run. These solutions were dominated by the exponential decay term of the GP model. The high likelihood of such solutions showed that the GP tends to favor models in which the data set for YZ Ceti is not dominated by the stellar rotation and in which there is no strong correlation between neighboring data points. The GP may also have a tendency to fit for such high α values due to the sampling of the RV data together with the intrinsic flexibility of the GP model. Nevertheless, we identified an interesting feature in Fig. 8a around a period of 68 days, with samples that have likelihood values similar to the samples in the plateau. For such periods, the GP allowed very low values for α_GP, which are more consistent with a rather stable quasi-periodic signal. We now tuned the quasi-periodic GP to model only the stellar rotational signal and sample this local maximum. In the following we show how this affects the log-evidence and derived planetary parameters, in particular the eccentricity of the planets.

The following two runs (b, c) show what happens if the GP is tuned to specifically model only the signal that can be directly attributed to the rotation period of the star. We constrained the prior of the α parameter to lower values, so that the GP will predominately fit more periodic signals. This tuning was physically motivated in our case, and forces the GP to not primarily model uncorrelated noise and act more like a sinusoid, but it still allows changes in the amplitude or phase shifts making it more flexible than a simple sinusoid. Since our prior for α excludes a number of high-likelihood samples distributed over the whole prior volume, it is expected that the log-evidence of a model with such a strong constraint is smaller than a model that can fit these solutions. If the GP is forced to predominantly fit the rotational signal of the star but no other unknown systematics (e.g., jitter), the median of the posterior of the derived eccentricity for the innermost planet and the upper boundary of the 95% density interval shift towards lower values (Table 6), which is exactly what we expected. However, there seems to be a “sweet spot” where further constraining of the α value leads to significantly lower maximum likelihood achieved within the sample distribution and to significantly worse log-evidence. This was the case when we further constrained α to values below 10⁻⁶ as we were then sampling only the tail of the contribution from the stellar rotation signal.

A widely used approach if a rotation period of the star is available from photometry is to use this information as an informative normal prior within the radial velocity fitting. We tested this approach in run d by adoptinga normal prior based on the photometric rotation period, and uncertainty derived from the R band as it has the largest uncertainty. In the specific case of YZ Ceti, this constraint has no influence on the derived planetary parameters compared to run a as the GP model still fits predominantly for high α values. The uncertainty of the GP rotational period equals the prior range for this hyperparameter, showcasing the same problem as for run a, namely that most of the posterior samples favor that the GP models to a lesser extent correlated effects and short-term noise, which seem to dominate over the contribution from the stellar rotational period. Only a constraint on α is able to change this behavior of the GP. It is reassuring that in the case of run c this constraint on α results in a distribution of the posterior of the rotation period for the GP that is close to a normal distribution and consistent with the photometric observations even without the need of an informative normal prior on the rotation period. We also tested whether any other period between 30 d and 100 d is consistent with such low α values by applying a broader uniform prior on the rotation period ranging from 30 d to 100 d as in the very first run that we performed, but keeping the constraint of run c on α. For this test we still found only one single mode for the posterior of the rotation period of the GP, which peaked at the same period as for run c, and the GP α-P diagram looked almost as in Fig. 8c. Additionally, the derived planetary parameters from run c were more consistent with a much simpler analysis that used sinusoids as activity models. Run c also led to more dynamically stable solutions than run a. For the physically motivated reasons mentioned before, we adopted run c as our final GP model for YZ Ceti.

Fig. 8 Gaussian-process alpha-period diagram (αGP vs. PGP) for four runs with different priors listed in Table 6. The color-coding shows a likelihood range of $\Delta \ln \mathcal{L}\le 10$ normalized to the highest achieved log-likelihood within all four runs, and can be compared between the different subplots and runs. Samples with a $\Delta \ln \mathcal{L}>10$ compared to the highest achieved likelihood are shown in gray.

Transit search with TESS

In addition to the extensive long-term photometry presented in Sect. 3.2, short-cadence observations from the TESS satellite (Ricker et al. 2015) were also available. YZ Ceti was observed in Sector 3 (Camera 1, CCD 1) from 20 September to 10 October, 2018. However, there were no TESS objects of interest (TOIs) listed on the TESS data alerts public website for this target. As in Luque et al. (2019) we performed an independent signal search applying the transit least-squares (TLS; Hippke & Heller 2019) algorithm on the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curve provided by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) on the Mikulski Archive for Space Telescopes (MAST)⁷. No signals were found in the process.

We therefore investigated whether we could rule out transits of the three known planets, and thus use this information to further constrain the minimum inclination of the system assuming coplanar orbits. YZ Ceti b has the highest transit probability (p ≈ R_*∕a_b) with p ≈ 4%, while YZ Ceti c and d have transit probabilities of p ≈ 3%. However, as it is also the smallest planet in the system, the transit signal of YZ Ceti b would be the most difficult to detect. To calculate the planetary radius we used the semi-empirical mass-radius relationship from Zeng et al. (2016). We derived R_b ≈ 0.93 R_⊕, R_c ≈ 1.05 R_⊕, and R_c ≈ 1.04 R_⊕. Assuming a circular orbit, the transit depth of YZ Ceti b would thus be about ${(r_b/R_{*})}^2\approx 0.29\%$ and the transit duration Δt ≈ 0.73 h. The standard deviation of the TESS PDCSAP light curve is 0.16%, which means that the planet would be easy to detect if it had a full transit. This was also confirmed by injecting fake box-transits into the data set and running the TLS signal search.

We concluded that the maximum inclination of the system must be such that a full transit of YZ Ceti b is excluded. This yields $i_{\text{max}}=\arccos \left(R_{*}/a_b\right)=87.43$ deg.

Fig. 9 Phase-folded RVs to the planetary periods based on the median posterior parameters of the stable solutions (Table 7).

6 N-body integrations

We tested the long-term stability of the YZ Ceti system by using the SyMBA N-body symplectic integrator (Duncan et al. 1998), which was modified to work in Jacobi coordinates (e.g., Lee & Peale 2003). Each posterior sample was integrated for a maximum of one million orbits of the inner planet with time steps of 0.02 d. However SyMBA reduces the time step during close encounters to ensure an accurate simulation. SyMBA also tests whether there are planet–planet or planet–star collisions, or planetary ejections, and if so interrupts the integration. A planet is considered lost and the system unstable if, at any time, (i) the mutual planet–planet separation is below the sum of their physical radii, assuming Jupiter mean density (i.e., planets undergo collision); (ii) the star–planet separation exceeds twice the initial semi-major axis of the outermost planet (r_max > 2a_d init), which we define as planetary ejection; (iii) the star–planet separation is below the physical stellar radius (R ≈ 0.00074 au), which we consider a collision with the star; and (iv) the semi-major axis receives a change of 30% compared to the initial value.These criteria efficiently detect unstable configurations and save computation time.

The inclination of YZ Ceti, i_max < 87.43 deg, was applied for the stability analysis under the assumption of co-planar orbits. In Table 7 we show the further constrained posterior parameters that we derived by allowing only solutions stable up to one million orbits of the inner planet. Compared to our previous estimates, we found lower values for the eccentricities. We also found that about 17.1% of the samples of our favored model (run c) were stable, showing that the compactness of the planetary orbits in the YZ Ceti system allows only a narrow range of planetary eccentricities. For comparison, the stable fraction of posterior samples without modeling the activity with a GP was only 1.7%. Figure 9 shows the phase plot of the RV models based on the posterior of the stable solutions. A corner plot of all the derived fit parameters for run c using juliet is displayed in Fig. A.1. In this plot we also highlight the stable sample (in blue). Additionally, we used ourstability analysis to search for a lower limit on the inclination. We found that for an inclination of i_min = 0.9 deg no sample solution was stable, providing a weak upper limit. Based on the stable solutions we derived some additional planetary parameters which are given in Table 8.

The median GP model of the stable solutions and its GLS periodogram are shown in Fig. 10. With the GLS periodogram we verified which periods are fitted by the GP. The plot also shows that the amplitude of the GP increases over the observational time span, indicating that the contribution of the radial velocity variations caused by the stellar rotation increased over time.

The dynamics of the YZ Ceti multiplanetary system has only been sparsely investigated since its discovery by Astudillo-Defru et al. (2017). From our 17.1% of the stable Kepler sample, we found that a fraction of about 22% of the solutions showed clear libration of the three-body Laplace angle Θ_L given as $${\Theta }_{\text{L}}=2{\lambda }_1-5{\lambda }_2+3{\lambda }_3,$$(3)

where λ₁, λ₂, and λ₃ are the mean longitudes of YZ Ceti b, c, and d. This result is in agreement with a purely theoretical result by Pichierri et al. (2019) based on the measured period ratios between the planets. However, Pichierri et al. (2019) used the period of 1.97 days for YZ Ceti b, which is not favored by our analysis.

A full dynamical analysis using self-consistent N-body fits to the RV data of this compact three-planet system is beyond the scope of this paper and will be carried out in a separate study (Stock et al., in prep.), for which the results of this paper will serve as a basis. However, we tested an alternative approach where possible unstable parameter combinations are penalized during the Kepler fit in order to push the fit towards the stable solutions.

Several possibilities for such on-the-fly tests are possible. Since higher eccentricities tend to disturb the planetary system, a smooth cutoff for eccentricities could be implemented. If guided by dynamical simulations, reasonable cutoff values can be obtained. The mutual distances of the planets in units of their Hill radii are an alternative. Again, minimum Hill radii separations could be inferred from dynamical calculations, but on a more general level; however, both approaches are ad hoc, and the choice of cutoff values will restrict the parameter distribution.

The third approach is therefore to use the angular momentum deficit (AMD) and the Hill stability formulated using the AMD to assess dynamical stability. In a series of papers, Laskar (1997, 2000), Laskar & Petit (2017), and Petit et al. (2018) developed an easy-to-use formulation for AMD stability. In short, AMD is the sum of planetary eccentricities and mutual inclinations weighted by planetary mass and orbital separation. Since this quantity is conserved it allows us to evaluate close possible encounters in the planetary system. In the last paper, Petit et al. (2018) derived a formulation of the Hill stability using AMD and compared its numerical N-body simulations. We implemented this AMD–Hill stability and used it to penalize unstable parameter combinations during an MCMC fitting procedure. As expected, the AMD–Hill stability criterion inhibits solutions with eccentricities that are toohigh. For YZ Ceti b and d, we find a flat distribution of eccentricities up to about 0.15 and 0.12, respectively,with a steep drop above these values. YZ Ceti c has a distribution that decreases continuously from 0 to 0.14. In Fig. 11 the eccentricities are selected from the full corner plot to show the posterior distributions for those parameters. Overall, the results from the AMD–Hill stability approach are consistent with our results from the stability analysis based on the posterior samples, indicating that the eccentricities of the planets must be lower than derived from a simple Keplerian three-planet fit to the data without modeling the activity.

Fig. 10 Median GP model and its GLS periodogram based on the stable posterior samples for YZ Ceti. In this plot the three-planet model is not included, and is subtracted from the RV data. We show only the densely sampled region between BJD 2 456 480 and BJD 2 458 520. The gray area indicates the interdecile range of the GP model.

Fig. 11 Selection from the corner plot showing the posterior distributions of the eccentricities from the AMD–Hill stability run. The blue line indicates the best fit. The vertical dashed lines denote the 68% posterior credibility intervals and the median.

7 Discussion and conclusions

Using additional 229 RV measurements of YZ Ceti compared to the discovery study (Astudillo-Defru et al. 2017), we constrained the true planetary configuration and resolved the aliases discussed in the literature (e.g., Robertson 2018). We achieved this by using the AliasFinder which uses a slightly modified version of the method suggested by Dawson & Fabrycky (2010) to disentangle aliases. The results from the AliasFinder analysis were in agreement with an analysis regarding the comparison of the maximum log-likelihood within the posterior samples of different three-planet model realizations. The most likely planetary configuration determined from our data and both analyses was a system of three planets at periods of 2.02 d, 3.06 d, and 4.66 d, which differs from previously published configurations for this system. The new configuration results in an almost optimal 3:2 commensurability of the periods of neighboring planets. We found no statistically significant evidence for an additional fourth companion orbiting YZ Ceti even though we analyzedmore than two times the number of RV measurements than did Astudillo-Defru et al. (2017). In particular, we found no sign of the tentative signal at 1.04 d, in contrast to the discovery study. However, we did observe variations in the RV data with a period around 68 d caused by the stellar rotation. In contrast to the discovery study, which adopted a stellar rotation period of 83 d, we found values of the rotation period of 68.4 ± 0.05 d and 68.5 ± 1.0 d based on combined V - and R-band photometric follow-up with a number of instruments, respectively.

YZ Ceti is an example of a relatively old star with a long rotational period quite distinct from the exoplanet periods, where the activity strongly influences the determination of the planetary parameters from the RV model. Due to precise photometry, we were able to link an apparent period in the RV residuals of a circular three-planet Keplerian fit to the rotational period of the star. After modeling the stellar rotational signal with a quasi-periodic Gaussian process, we derived a lower eccentricity for the innermost planet than without modeling activity at all. This result is also more consistent with stability constraints that apply to this compact system. We found very good agreement between the photometric rotation period and the rotation period derived by the GP solely from the RV data. We observed only a small region where the quasi-periodic GP allowed low values for the inverse-lengthscale α, consistent with a rather stable periodic signal. This small region was consistent with the estimates of the stellar rotation period from photometry.

Interestingly, the second harmonic of the orbital period of YZ Ceti b, which is related to the eccentricity, is very close to an alias of the rotation period. Thus, incorrectly modeled stellar activity RV modulations can cause deviations from a sinusoid in the reflex RV curve of YZ Ceti b. The Keplerian fit to YZ Ceti b accommodates this by fitting an eccentric orbit. This may explain the surprisingly strong influence of the modeling of the rotational variations on the derived eccentricity of YZ Ceti b.

We searched for transits using TESS light curves. From their non-detection we derived an upper limit to the inclination of the system of i_max = 87.43 deg. Applying the criterion of long-term stability, we were able to reduce the uncertainties of the planetary parameters. We also determined a weak lower limit for the inclination of the planets, which is i_min = 0.9 deg. Additionally, we noted that for 22% of the stable orbital integrations the three-body resonance angle librates, so it is possible that a resonant chain was established during the formation of the ultra-compact YZ Ceti system.

Overall, the detailed analysis outlined within this work shows how different novel techniques can help to constrain the architecture of systems hosted by active stars.

Acknowledgements

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 de Calar Alto (CAHA, Almería, Spain). CARMENES is funded by the German Max-Planck-Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF FICTS – 2011 – 02 funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Insitut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano - Alemán), with additional contributions by the Spanish Ministry of Economy, the German Science Foundation through the Major Research Instrumentation Programme 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. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 37 / 935 – 1 FUGG. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia, Innovacióny Universidades and the European FEDER/ERF funds through projects AYA2015-69350-C3-2-P, AYA2016-79425-C3-1/2/3-P, ESP2017-87676-C5-2-R, ESP2017-87143-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ísicade Andalucía (SEV-2017-0709), and Centro de Astrobiología (MDM-2017-0737), and the Generalitat de Catalunya/CERCA programme. This work is supported by the Science and Technology Facilities Council [ST/M001008/1 and ST/P000584/1]. J.S.J. acknowledges support by Fondecyt grant 1161218 and partial support by CATA-Basal (PB06, CONICYT). Z.M.B acknowledges funds from CONICYT-FONDECYT/Chile Postdoctorado 3180405. M.H.L. was supported in part by Hong Kong RGC grant HKU 17305618. T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. Data were partly collected with the robotic 40-cm telescope ASH2 at the SPACEOBS observatory (San Pedro de Atacama, Chile) and the 90-cm telescope at Sierra Nevada Observatory (Granada, Spain) both operated by the Instituto de Astrofífica de Andalucía (IAA). This work has made use of the Exo-Striker tool (Trifonov 2019).

Appendix A Cornerplot of stable solutions

Fig. A.1 Corner plot of the planetary parameters. The gray areas indicate the different sigma levels of the juliet samples with the GP model (run c in Table 6). The blue points show the distribution of samples stable over 1 million orbits of the innermost planet (approximately 5557 yr). Error bars denote the 68% posterior credibility intervals.

Appendix B RV data

Table B.1 is available at the CDS.

Appendix C Additional table

References

All Tables

All Figures

Fig. 1 Top three rows: nightly binned photometric time series (top: all data in the V band, middle: zoom to ASH2 and SNO data in the V band, bottom: all data in the R band) and aphase plot to the determined rotation period. Bottom row: periodograms of the different instruments in the V band (left), R band (middle), and the combination of all instruments in each band (right). The periodograms for the analysis on individual instruments are color-coded (blue: ASH2, brown: SNO, red: ASAS, green: ASAS-SN, dark green: MEarth T11, purple:MEarth T12), while the combined periodograms are plotted in black. The solid line represents the combined V -band periodogramand the dashed line the combined R-band periodogram. For the combined periodograms we show the FAP level of 0.001 (green solid line for V band and red dashed line for R band). The vertical black line in each periodogram represents the determined rotation period of 68.4 d and 68.5 d, respectively.

In the text

	Fig. 2 GLS periodograms of several activity indicators in the CARMENES and HARPS data. The periodograms are separated and show different frequency regimes to better display the significant peaks within the low-frequency regime. The red dashed line shows the photometric rotation period and the black solid lines highlight the planetary frequencies.
In the text

	Fig. 3 Radial velocity data and final stable fit including the Gaussian process model for the activity signal (see Sects. 5 and 6).
In the text

	Fig. 4 Generalized Lomb-Scargle periodograms of the data and different fit models from Table C.1, and the final GLS periodogram after fitting for the activity with a GP (see Sect. 5).
In the text

	Fig. 5 S-BGLS periodograms after subtracting the three planetary signals. Left: around the frequency of the tentative signal mentioned by Astudillo-Defru et al. (2017), middle: 29.36-day signal visible as remnant power in a circular three-planetmodel, right: 68-day signal attributed to the stellar rotation.
In the text

	Fig. 6 Periodogram of the spectral window function for the data used in our analysis.
In the text

Fig. 7 Alias tests for the periods of 3.06 d (a), 4.66 d (b), and 2.02 d (c). In each plot,each row corresponds to one set of simulations. The frequency of the injected signal is indicated by a vertical blue dashedline. From 1000 simulated data sets each, the median of the obtained periodograms (black solid line), the interquartile range, and the ranges of 90 and 99% (shades of gray) are shown. For comparison, the periodogram of the observed data is plotted as a red solid line. Additionally, the angular mean of the phase of each peak and its standard deviation of the simulated periodograms are shown as clock diagrams (black line and grays) and can be compared to the phase of thepeak in the observed periodogram (red line).

In the text

Fig. 8 Gaussian-process alpha-period diagram (αGP vs. PGP) for four runs with different priors listed in Table 6. The color-coding shows a likelihood range of $\Delta \ln \mathcal{L}\le 10$ normalized to the highest achieved log-likelihood within all four runs, and can be compared between the different subplots and runs. Samples with a $\Delta \ln \mathcal{L}>10$ compared to the highest achieved likelihood are shown in gray.

In the text

	Fig. 9 Phase-folded RVs to the planetary periods based on the median posterior parameters of the stable solutions (Table 7).
In the text

	Fig. 10 Median GP model and its GLS periodogram based on the stable posterior samples for YZ Ceti. In this plot the three-planet model is not included, and is subtracted from the RV data. We show only the densely sampled region between BJD 2 456 480 and BJD 2 458 520. The gray area indicates the interdecile range of the GP model.
In the text

	Fig. 11 Selection from the corner plot showing the posterior distributions of the eccentricities from the AMD–Hill stability run. The blue line indicates the best fit. The vertical dashed lines denote the 68% posterior credibility intervals and the median.
In the text

	Fig. A.1 Corner plot of the planetary parameters. The gray areas indicate the different sigma levels of the juliet samples with the GP model (run c in Table 6). The blue points show the distribution of samples stable over 1 million orbits of the innermost planet (approximately 5557 yr). Error bars denote the 68% posterior credibility intervals.
In the text