Issue 
A&A
Volume 636, April 2020



Article Number  A119  
Number of page(s)  19  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201936732  
Published online  01 May 2020 
The CARMENES search for exoplanets around M dwarfs
Characterization of the nearby ultracompact multiplanetary system YZ Ceti^{★}
^{1}
Landessternwarte, Zentrum für Astronomie der Universität Heidelberg,
Königstuhl 12,
69117
Heidelberg,
Germany
email: sstock@lsw.uniheidelberg.de
^{2}
MaxPlanckInstitut für Astronomie,
Königstuhl 17,
69117
Heidelberg, Germany
^{3}
Institut für Astrophysik, GeorgAugustUniversität,
FriedrichHundPlatz 1,
37077
Göttingen, Germany
^{4}
Centro de Astrobiología (CSICINTA), ESAC,
Camino bajo del castillo s/n,
28692
Villanueva de la Cañada,
Madrid, Spain
^{5}
School of Physics and Astronomy, Queen Mary University of London,
327 Mile End Road,
London,
E1 4NS, UK
^{6}
Instituto de Astrofísica de Andalucía (IAACSIC),
Glorieta de la Astronomía s/n,
18008
Granada, Spain
^{7}
Institut de Ciències de l’Espai (ICE, CSIC),
Campus UAB, C/Can Magrans s/n,
08193
Bellaterra, Spain
^{8}
Institut d’Estudis Espacials de Catalunya (IEEC),
08034
Barcelona, Spain
^{9}
School of Physical Sciences, The Open University,
Milton Keynes,
MK7 6AA, UK
^{10}
Departamento de Astronomía, Universidad de Chile,
Camino El Observatorio 1515,
Las Condes,
Casilla 36D,
Santiago, Chile
^{11}
Instituto de Astrofísica de Canarias (IAC),
38205
La Laguna,
Tenerife, Spain
^{12}
Departamento de Astrofísica, Universidad de La Laguna (ULL),
38206
La Laguna,
Tenerife, Spain
^{13}
Physikalisches Institut, Universitaet Bern,
Gesellschaftsstrasse 6,
3012
Bern, Switzerland
^{14}
Departamento de Física de la Tierra y Astrofísica & IPARCOSUCM (Instituto de Física de Partículas y del Cosmos de la UCM), Facultad de Ciencias Físicas, Universidad Complutense de Madrid,
28040
Madrid, Spain
^{15}
Department of Exploitation and Exploration of Mines, University of Oviedo,
Oviedo, Spain
^{16}
Thüringer Landessternwarte Tautenburg,
Sternwarte 5,
07778
Tautenburg, Germany
^{17}
Departamento de Astronomia, Universidad de Chile,
Casilla 36D,
Santiago, Chile
^{18}
Centro de Astrofísica y Tecnologías Afines (CATA),
Casilla 36D,
Santiago, Chile
^{19}
Centre for Astrophysics Research, University of Hertfordshire,
Hatfield
AL10 9AB, UK
^{20}
Department of Earth Sciences and Department of Physics, The University of Hong Kong,
Pokfulam Road,
Hong Kong
^{21}
Observatorio de Calar Alto, Sierra de los Filabres,
04550
Gérgal,
Almería, Spain
Received:
18
September
2019
Accepted:
20
January
2020
Context. The nearby ultracompact multiplanetary system YZ Ceti consists of at least three planets, and a fourth tentative signal. The orbital period of each planet is the subject of discussion in the literature due to strong aliasing in the radial velocity data. The stellar activity of this M dwarf also hampers significantly the derivation of the planetary parameters.
Aims. With an additional 229 radial velocity measurements obtained since the discovery publication, we reanalyze the YZ Ceti system and resolve the alias issues.
Methods. We use model comparison in the framework of Bayesian statistics and periodogram simulations based on a method by Dawson and Fabrycky to resolve the aliases. We discuss additional signals in the RV data, and derive the planetary parameters by simultaneously modeling the stellar activity with a Gaussian process regression model. To constrain the planetary parameters further we apply a stability analysis on our ensemble of Keplerian fits.
Results. We find no evidence for a fourth possible companion. We resolve the aliases: the three planets orbit the star with periods of 2.02 d, 3.06 d, and 4.66 d. We also investigate an effect of the stellar rotational signal on the derivation of the planetary parameters, in particular the eccentricity of the innermost planet. Using photometry we determine the stellar rotational period to be close to 68 d and we also detect this signal in the residuals of a threeplanet fit to the RV data and the spectral activity indicators. From our stability analysis we derive a lower limit on the inclination of the system with the assumption of coplanar orbits which is i_{min} = 0.9 deg. From the absence of a transit event with TESS, we derive an upper limit of the inclination of i_{max} = 87.43 deg.
Conclusions. YZ Ceti is a prime example of a system where strong aliasing hindered the determination of the orbital periods of exoplanets. Additionally, stellar activity influences the derivation of planetary parameters and modeling them correctly is important for the reliable estimation of the orbital parameters in this specific compact system. Stability considerations then allow additional constraints to be placed on the planetary parameters.
Key words: techniques: radial velocities / stars: individual: YZ Ceti / stars: latetype / planets and satellites: terrestrial planets / planets and satellites: dynamical evolution and stability / planetary systems
Table B.1 is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/636/A119.
© 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 spectroscopydetected systems had stellar masses lower than 0.3 M_{⊙}. Such a lowmass star is YZ Ceti (GJ 54.1), which was reported to host three planets (AstudilloDefru et al. 2017). It is the closest multiplanetary system to our Solar System published so far. AstudilloDefru et al. (2017) announced that YZ Ceti is orbited by at least three Earthmass 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^{−1}) 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 (AstudilloDefru 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 AstudilloDefru et al. (2017). In particular, the signal of planet c (P = 3.06 d) had a strong alias at 0.75 days.
AstudilloDefru et al. (2017) also mentioned a fourth tentative signal slightly above a oneday 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 ultracompact planetary systems around stars similar to YZ Ceti should be common, where the planets are typically locked in long resonant chains, exhibiting both twobody and threebody 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 longterm 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^{1} 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 nearinfrared (NIR) the 0.96–1.71 μm range with spectral resolution of R = 81 800 (Quirrenbach et al. 2014). We obtained 111 highresolution 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 FabryPé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 signaltonoise 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 leastsquares fitting. The radial velocities were corrected for barycentric motion, secular perspective acceleration, instrumental drift, and nightly zero points (see Trifonov et al. 2018; TalOr 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^{−1} and a root mean square (rms) of 3.00 m s^{−1} 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^{−1} and an rms of 7.17 m s^{−1} around the mean value. Due to the small amplitudes of the planetary signals of less than 2 m s^{−1} 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 subm s^{−1} precision. We retrieved 334 highresolution 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^{−1}. 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 postfiber data and fitted an offset for it. For the prefiber upgrade data sets we achieved a median internal uncertainty of 1.92 m s^{−1} and an rms of 2.88 m s^{−1} around the mean value. For the postfiber upgrade data sets we achieved a median internal uncertainty of 1.86 m s^{−1} and an rms of 3.72 m s^{−1} around the mean value.
2.3 Photometry
Details of observations from five photometric facilties are given below.
ASAS
The AllSky Automated Survey (Pojmański 1997) has been monitoring the entire southern and part of the northern sky since 2000. We retrieved 461 ASAS3 measurements of YZ Ceti taken between November 2000 and November 2009.
ASASSN
The AllSky Automated Survey for Supernovae (ASASSN) (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^{2} at different observatories worldwide. We extracted about six years of photometric observations (2013–2019) in the V band from the ASASSN archive^{2} for YZ Ceti.
MEarth
The MEarth project (Berta et al. 2012) is an allsky 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 InterAmerican 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^{2}. MEarth generally uses an RG715^{3} longpass filter, except for the 20102011 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^{4}). 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^{2}. 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 RitcheyChré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^{5}. 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 StefanBoltzmann law. With a linear massradius 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ésContreras (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ésContreras 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 Xray emission, which is consistent with the star kinematics.
3.2 Photometric analysis
AstudilloDefru et al. (2017) estimated the rotation period of YZ Ceti to be 83 d by analyzing ASAS photometry and the FWHM of the crosscorrelation 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 ASASSN 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 AstudilloDefru 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, ASASSN, 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 LombScargle (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 χ^{2} 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 Rband 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 falsealarm 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 AstudilloDefru et al. (2017). Fitting a sinusoid for either the 68day signal or the 80day signal removes the other, a strong sign of aliasing. In all data sets the 80day signal can be reproduced by strong aliases due to annual sampling effects in the window function together with the 68day periodicity. There are no strong signs of the 80day periodicity in the combined Rband 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.
Stellar parameters of YZ Ceti.
Highest GLS peak (P_{high}) and alternative GLS peak consistent with a common rotation period (P_{rot}) ^{(a)}.
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^{−1}, 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^{−1}, well below the directly estimated upper limit of 2 km s^{−1}. 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 crosscorrelation 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 HARPSPOST 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 AstudilloDefru 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 (ForemanMackey 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 logevidence of the corresponding models. The resulting logevidence 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 logevidence differs by . 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 [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 AstudilloDefru 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 Nbody integrator (see Sect. 5 and Table 7). Therefore, we compared the Bayesian logevidence 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 () or fixes eccentricity to zero for all planets (). As a result, there is moderate evidence to fit eccentric orbits for planets b and d (Trotta 2008). The residual periodogram of a threeplanet 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 stackedBayesian GLS periodogram (sBGLS) 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 sBGLS periodograms to their respective minimum values.
We found that the 28.37day signal (f ≈ 0.035 d^{−1}) 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 68day 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 threeplanet circular fit. After fitting a threeplanet 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^{−1}, 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 AstudilloDefru et al. (2017) at 1.04 d (f ≈ 0.962 d^{−1}) has decreased in significance (see also Fig. 5, left). In contrast to these signals, our sBGLS 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 colorcoded (blue: ASH2, brown: SNO, red: ASAS, green: ASASSN, 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 Rband 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 lowfrequency 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. AstudilloDefru 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 loglikelihood, similar to Robertson (2018), of different realizations of periods for the threeplanet 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 instead of . The loglikelihood is more robust to changes in the prior than the logevidence when comparing the same model with equal number of free parameters but different realizations. The achieved maximum loglikelihoods 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 (, 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 . 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 firstorder aliases. The aliases were thereby calculated using the equation (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 firstorder daily aliases of f_{alias} = f_{p} − f_{s} and f_{alias} = f_{p} + f_{s} respectively (f_{s} = 1.0027 d^{−1}). 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 aliastesting on github^{6} (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^{−1}) 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.06day signal is injected into our simulation compared to when we use its aliases. Injecting the alias at 1.48 d (f ≈ 0.676 d^{−1}, 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^{−1}). Injecting the 0.75day 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^{−1}) fits the peak heights and phases of the data periodogram well. Injecting the signal of 1.27 d (f ≈ 0.788 d^{−1}; see Fig. 7b, middle panel) did not reproduce the peak at 0.82 days (f ≈ 1.218 d^{−1}) and resultedin larger differences in the phase values of the other two peaks. Instead, the 0.82day signal (Fig. 7b, lower panel) did not reproduce the 1.27day peak. We concluded that YZ Ceti d orbits at 4.66 d as stated by AstudilloDefru 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^{−1}) published by AstudilloDefru et al. (2017) performs worse than the alternative period of 2.02 d (f ≈ 0.495 d^{−1}; see Fig. 7c, middle panel). The 1.97day signal did not reproduce the peak at 0.67 d (f ≈ 1.497 d^{−1}) 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.67daysignal (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 AstudilloDefru 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 threeplanet 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). 
Bayesian logevidence for models of different number of planets and their logevidence difference to the best model.
Fig. 4 Generalized LombScargle 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 SBGLS periodograms after subtracting the three planetary signals. Left: around the frequency of the tentative signal mentioned by AstudilloDefru et al. (2017), middle: 29.36day signal visible as remnant power in a circular threeplanetmodel, right: 68day signal attributed to the stellar rotation. 
Maximum achieved loglikelihood for different threeplanet models and their difference to the best model.
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). 
Planetary orbital frequencies f_{p} and their firstorder aliases for a sampling frequency of f_{s} = 1.0027d^{−1}.
5 Simultaneous fitting of stellar activity and Keplerians
As in AstudilloDefru 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 expsinsquared kernel multiplied with a squaredexponential kernel (Ambikasaran et al. 2015). This kernel has the form (2)
where is the amplitude of the GP component given in m s^{−1}, α_{i} is the inverse length scale of the GP exponential component given in d^{−2}, Γ_{i} is the amplitude of GP sinesquared component given in m s^{−1}, P_{rot,i} the period of the GP quasiperiodic component given in d, and τ is the timelag. The α value is a measure of the strength of the exponential decay of the quasiperiodic kernel. A lower α describes a more stable periodic signal and stronger correlation between data points. The quasiperiodic 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 quasiperiodic kernel and the celerite Simple Harmonic Oscillator (SHO) kernel (ForemanMackey et al. 2017). Both of these had significantly worse evidence than the quasiperiodic 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 highdimensional 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 . 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^{0} d^{−2} and 10^{−4} d^{−2}. 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 loglikelihood, and the eccentricity of planet b as this parameter is rather sensitive to the modeling of the stellar activity. The Bayesian logevidence of all these runs was significantly better than a simple threeplanet Keplerian fit to the data. In Fig. 8 we show the scatter plots of the sampled α values of the quasiperiodic 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 logevidence while increasing the number of parameters. The derived planetary parameters and remaining nonamplitude 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 logevidence, 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^{−2} to 1 d^{−2}, corresponding to decay timescales of several days (). The plateau spanned the range between 10^{0} d^{−2} and 10^{−4} d^{−2} 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 quasiperiodic signal. We now tuned the quasiperiodic GP to model only the stellar rotational signal and sample this local maximum. In the following we show how this affects the logevidence 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 highlikelihood samples distributed over the whole prior volume, it is expected that the logevidence 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 logevidence. This was the case when we further constrained α to values below 10^{−6} 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 shortterm 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.
Four runs with different GP priors done with juliet for modeling the activity of YZ Ceti.
Fig. 8 Gaussianprocess alphaperiod diagram (α_{GP} vs. P_{GP}) for four runs with different priors listed in Table 6. The colorcoding shows a likelihood range of normalized to the highest achieved loglikelihood within all four runs, and can be compared between the different subplots and runs. Samples with a compared to the highest achieved likelihood are shown in gray. 
Posterior parameters of fits obtained for YZ Ceti using juliet.
Transit search with TESS
In addition to the extensive longterm photometry presented in Sect. 3.2, shortcadence 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 leastsquares (TLS; Hippke & Heller 2019) algorithm on the Presearch 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)^{7}. 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 semiempirical massradius 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 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 boxtransits 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 deg.
Fig. 9 Phasefolded RVs to the planetary periods based on the median posterior parameters of the stable solutions (Table 7). 
6 Nbody integrations
We tested the longterm stability of the YZ Ceti system by using the SyMBA Nbody 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 semimajor 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 semimajor 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 coplanar 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 AstudilloDefru 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 threebody Laplace angle Θ_{L} given as (3)
where λ_{1}, λ_{2}, and λ_{3} 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 selfconsistent Nbody fits to the RV data of this compact threeplanet 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 onthefly 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 easytouse 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 Nbody 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 threeplanet fit to the data without modeling the activity.
Derived planetary parameters obtained for YZ Ceti b, c, and d.
Fig. 10 Median GP model and its GLS periodogram based on the stable posterior samples for YZ Ceti. In this plot the threeplanet 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 (AstudilloDefru 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 loglikelihood within the posterior samples of different threeplanet 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 AstudilloDefru 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 Rband photometric followup 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 threeplanet Keplerian fit to the rotational period of the star. After modeling the stellar rotational signal with a quasiperiodic 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 quasiperiodic GP allowed low values for the inverselengthscale α, 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 nondetection we derived an upper limit to the inclination of the system of i_{max} = 87.43 deg. Applying the criterion of longterm 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 threebody resonance angle librates, so it is possible that a resonant chain was established during the formation of the ultracompact 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 / 41. CARMENES is an instrument for the Centro Astronómico HispanoAlemán de Calar Alto (CAHA, Almería, Spain). CARMENES is funded by the German MaxPlanckGesellschaft (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 (MaxPlanckInstitut 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 BadenWü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 BadenWü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 AYA201569350C32P, AYA201679425C31/2/3P, ESP201787676C52R, ESP201787143R and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Canarias (SEV20150548), Instituto de Astrofísicade Andalucía (SEV20170709), and Centro de Astrobiología (MDM20170737), 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 CATABasal (PB06, CONICYT). Z.M.B acknowledges funds from CONICYTFONDECYT/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 40cm telescope ASH2 at the SPACEOBS observatory (San Pedro de Atacama, Chile) and the 90cm 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 ExoStriker 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
Priors used within juliet to model the YZ Ceti multiplanetary system.
References
 AlonsoFloriano, F. J., Morales, J. C., Caballero, J. A., et al. 2015, A&A, 577, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252 [NASA ADS] [CrossRef] [Google Scholar]
 AstudilloDefru, N., Díaz, R. F., Bonfils, X., et al. 2017, A&A, 605, L11 [Google Scholar]
 Berta, Z. K., Irwin, J., Charbonneau, D., Burke, C. J., & Falco, E. E. 2012, AJ, 144, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Caballero, J. A., CortésContreras, M., AlonsoFloriano, F. J., et al. 2016a, 19th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 148 [Google Scholar]
 Caballero, J. A., Guàrdia, J., López del Fresno, M., et al. 2016b, Proc. SPIE, 9910, 99100E [Google Scholar]
 Coleman, G. A. L., Leleu, A., Alibert, Y., & Benz, W. 2019, A&A, 631, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CortésContreras, M. 2016, PhD thesis, Universidad Complutense de Madrid, Madrid, Spain [Google Scholar]
 Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Díez Alonso, E., Caballero, J. A., Montes, D., et al. 2019, A&A, 621, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dreizler, S., Jeffers, S. V., Rodríguez, E., et al. 2020, MNRAS, 493, 536 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Engle, S. G., & Guinan, E. F. 2017, ATel, 10678 [Google Scholar]
 Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Agol, E., Angus, R., & Ambikasaran, S. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gliese, W., & Jahreiß, H. 1995, VizieR Online Data Catalog: V/70A [Google Scholar]
 Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Husser, T.O., Wendevon Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2017, ATel, 10643 [Google Scholar]
 Jeffreys, H. 1946, Proc. R. Soc. London Ser. A, 186, 453 [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E [Google Scholar]
 Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502 [NASA ADS] [CrossRef] [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [NASA ADS] [CrossRef] [Google Scholar]
 Kunkel, W. E. 1970, IBVS, 442, 1 [NASA ADS] [Google Scholar]
 Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J. 1997, A&A, 317, L75 [NASA ADS] [Google Scholar]
 Laskar, J. 2000, Phys. Rev. Lett., 84, 3240 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Laskar, J., & Petit, A. C. 2017, A&A, 605, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Lo Curto, G., Pepe, F., Avila, G., et al. 2015, The Messenger, 162, 9 [NASA ADS] [Google Scholar]
 Luque, R., Pallé, E., Kossakowski, D., et al. 2019, A&A, 628, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Montes, D., LópezSantiago, J., Gálvez, M. C., et al. 2001, MNRAS, 328, 45 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mortier, A., & Collier Cameron, A. 2017, A&A, 601, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mortier, A., Faria, J. P., Correia, C. M., Santerne, A., & Santos, N. C. 2015, A&A, 573, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Passegger, V. M., Reiners, A., Jeffers, S. V., et al. 2018, A&A, 615, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petit, A. C., Laskar, J., & Boué, G. 2018, A&A, 617, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pichierri, G., Batygin, K., & Morbidelli, A. r. 2019, A&A, 625, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pojmański, G. 1997, Acta Astron., 47, 467 [NASA ADS] [Google Scholar]
 Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, Proc. SPIE, 9147, 91471F [Google Scholar]
 Reiners, A., Ribas, I., Zechmeister, M., et al. 2018, A&A, 609, L5 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum., and Systems, 1, 014003 [Google Scholar]
 Roberts, D. H., Lehar, J., & Dreher, J. W. 1987, AJ, 93, 968 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, P. 2018, ApJ, 864, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Rodríguez, E., García, J. M., Costa, V., et al. 2010, MNRAS, 408, 2149 [NASA ADS] [CrossRef] [Google Scholar]
 Schweitzer, A., Passegger, V. M., Cifuentes, C., et al. 2019, A&A, 625, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakhovskaya,N. I. 1995, VizieR Online Data Catalog: II/055 [Google Scholar]
 Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Speagle, J. S. 2020, MNRAS, 493, 3132 [NASA ADS] [CrossRef] [Google Scholar]
 Stock, S., & Kemmer, J. 2020, J. Open Source Softw. 5, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Suárez Mascare no, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 TalOr, L., Trifonov, T., Zucker, S., Mazeh, T., & Zechmeister, M. 2019, MNRAS, 484, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Trifonov, T. 2019, Astrophysics Source Code Library [record ascl:1906.004] [Google Scholar]
 Trifonov, T., Kürster, M., Zechmeister, M., et al. 2018, A&A, 609, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trotta, R. 2008, Contemp. Phys., 49, 71 [Google Scholar]
 Tuomi, M., Jones, H. R. A., AngladaEscudé, G., et al. 2019, AAS J., submitted [arXiv:1906.04644] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12 [CrossRef] [Google Scholar]
 Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Highest GLS peak (P_{high}) and alternative GLS peak consistent with a common rotation period (P_{rot}) ^{(a)}.
Bayesian logevidence for models of different number of planets and their logevidence difference to the best model.
Maximum achieved loglikelihood for different threeplanet models and their difference to the best model.
Planetary orbital frequencies f_{p} and their firstorder aliases for a sampling frequency of f_{s} = 1.0027d^{−1}.
Four runs with different GP priors done with juliet for modeling the activity of YZ Ceti.
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 colorcoded (blue: ASH2, brown: SNO, red: ASAS, green: ASASSN, 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 Rband 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 lowfrequency 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 LombScargle 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 SBGLS periodograms after subtracting the three planetary signals. Left: around the frequency of the tentative signal mentioned by AstudilloDefru et al. (2017), middle: 29.36day signal visible as remnant power in a circular threeplanetmodel, right: 68day 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 Gaussianprocess alphaperiod diagram (α_{GP} vs. P_{GP}) for four runs with different priors listed in Table 6. The colorcoding shows a likelihood range of normalized to the highest achieved loglikelihood within all four runs, and can be compared between the different subplots and runs. Samples with a compared to the highest achieved likelihood are shown in gray. 

In the text 
Fig. 9 Phasefolded 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 threeplanet 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.