Issue 
A&A
Volume 636, April 2020



Article Number  L6  
Number of page(s)  17  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/201937254  
Published online  16 April 2020 
Letter to the Editor
The SOPHIE search for northern extrasolar planets
XVI. HD 158259: A compact planetary system in a near3:2 mean motion resonance chain
^{1}
Observatoire Astronomique de l’Université de Genéve, 51 Chemin des Maillettes, 1290 Versoix, Switzerland
email: nathan.hara@unige.ch
^{2}
Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
^{3}
Center of Excellence in Information Systems, Tennessee State University, Nashville, TN 37209, USA
^{4}
Observatoire de HauteProvence, CNRS, Aix Marseille Universitté, Institut Pythtéas UMS 3470, 04870 SaintMichell’Observatoire, France
^{5}
Departamento de Matemática y Física Aplicadas, Universidad Católica de la Santísima Concepción, Alonso de Rivera 2850, Concepción, Chile
^{6}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
^{7}
Institut d’Astrophysique de Paris, UMR7095 CNRS, Universitté Pierre & Marie Curie, 98bis boulevard Arago, 75014 Paris, France
^{8}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150762 Porto, Portugal
^{9}
Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Buenos Aires, Argentina
^{10}
CONICET – Universidad de Buenos Aires, Instituto de Astronomía y Física del Espacio (IAFE), Buenos Aires, Argentina
^{11}
Instituto de Astrofésica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, Macul, Santiago, Chile
^{12}
Millennium Institute for Astrophysics, Macul, Chile
^{13}
CanadaFranceHawaii Telescope Corporation, 651238 Mamalahoa Hwy, Kamuela, HI 96743, USA
^{14}
Departamento de Fésica e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169007 Porto, Portugal
^{15}
Department of Physics, University of Warwick, Coventry CV4 7AL, UK
^{16}
Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, UK
Received:
4
December
2019
Accepted:
10
March
2020
Aims. Since 2011, the SOPHIE spectrograph has been used to search for Neptunes and superEarths in the northern hemisphere. As part of this observational program, 290 radial velocity measurements of the 6.4 V magnitude star HD 158259 were obtained. Additionally, TESS photometric measurements of this target are available. We present an analysis of the SOPHIE data and compare our results with the output of the TESS pipeline.
Methods. The radial velocity data, ancillary spectroscopic indices, and groundbased photometric measurements were analyzed with classical and ℓ_{1} periodograms. The stellar activity was modeled as a correlated Gaussian noise and its impact on the planet detection was measured with a new technique.
Results. The SOPHIE data support the detection of five planets, each with m sin i ≈ 6 M_{⊕}, orbiting HD 158259 in 3.4, 5.2, 7.9, 12, and 17.4 days. Though a planetary origin is strongly favored, the 17.4 d signal is classified as a planet candidate due to a slightly lower statistical significance and to its proximity to the expected stellar rotation period. The data also present low frequency variations, most likely originating from a magnetic cycle and instrument systematics. Furthermore, the TESS pipeline reports a significant signal at 2.17 days corresponding to a planet of radius ≈1.2 R_{⊕}. A compatible signal is seen in the radial velocities, which confirms the detection of an additional planet and yields a ≈2 M_{⊕} mass estimate.
Conclusions. We find a system of five planets and a strong candidate near a 3:2 mean motion resonance chain orbiting HD 158259. The planets are found to be outside of the two and three body resonances.
Key words: planets and satellites: detection / planets and satellites: dynamical evolution and stability / planets and satellites: fundamental parameters / planets and satellites: formation / methods: statistical / techniques: radial velocities
© ESO 2020
1. Introduction
Transit surveys have unveiled several multiplanetary systems where the planets are tightly spaced and close to low order mean motion resonances (MMRs). For instance, Kepler80 (Xie 2013; Lissauer et al. 2014; Shallue & Vanderburg 2018), Kepler223 (Borucki et al. 2011; Mills et al. 2016), and TRAPPIST1 (Gillon et al. 2016; Luger et al. 2017) present 5, 4, and 7 planets, respectively, in such configurations. These systems are often qualified as compact, in the sense that any two subsequent planets have a period ratio below 2. Compact, near resonant configurations could be the result of a formation scenario where the planets encounter dissipation in the gas disk, are locked in resonance, and then migrate inwards before potentially leaving the resonance (e.g., Terquem & Papaloizou 2007; MacDonald et al. 2016; Izidoro et al. 2017).
Near resonant, compact systems are detectable by radial velocity (RV), as demonstrated by follow up observations of transits (Lopez et al. 2019). However, such detections with only RV are rare: HD 40307 (Mayor et al. 2009) and HD 215152 (Delisle et al. 2018) both have three planets near 2:1 – 2:1 and 5:3 – 3:2 configurations, respectively.
In the present work, we analyze the 290 SOPHIE radial velocity measurements of HD 158259. We detect several signals, which are compatible with a chain of near resonant planets. The signals have an amplitude in the 1 − 3 m s^{−1} range. At this level, in order to confirm their planetary origin, it is critical to consider whether these signals could be due to the star or to instrument systematics. To this end, we include the following data sets in our analysis: the bisector span and derived from the spectra as well as groundbased photometric data. The periodicity search is performed with a ℓ_{1} periodogram (Hara et al. 2017) including a correlated noise model, selected with new techniques. The results are compared to those of a classical periodogram (Baluev 2008).
Furthermore, HD 158259 has been observed in sector 17 of the TESS mission (Ricker et al. 2014). The results of the TESS reduction pipeline (Jenkins et al. 2010, 2016) are included in our analysis.
The data support the detection of six planets close to a 3:2 MMR chain, with a lower detection confidence for the outermost one. The orbital stability of the resulting system is checked with numerical simulations, and we discuss whether the system is in or out of the two and threebody resonances.
The Letter is structured as follows. The data and its analysis are presented in Sects. 2 and 3, respectively. The study of the system dynamics is presented in Sect. 4, and we conclude in Sect. 5.
2. Data
2.1. HD 158259
HD 158259 is a G0 type star in the northern hemisphere with a V magnitude of 6.4. The known stellar parameters are reported in Table 1. The stellar rotation period is not known precisely, but it can be estimated. The median , which was obtained from SOPHIE measurements, is −4.8 ± 0.1. With the empirical relationship of Mamajek & Hillenbrand (2008), this translates to an estimated rotation period of 18 ± 5 days. Additionally, the SOPHIE RV data give v sin i = 2.9 ± 1 km s^{−1} (see Boisse et al. 2010). Assuming i = 90° and taking the Gaia radius estimate of 1.21 R_{⊙}, the v sin i estimation yields a rotation period of ≈20 ± 7 days.
Known stellar parameters of HD 158259.
2.2. SOPHIE radial velocities
SOPHIE is an échelle spectrograph mounted on the 193cm telescope of the HauteProvence Observatory (Bouchy et al. 2011). Several surveys have been conducted with SOPHIE, including a moderate precision survey (3.5–7 m s^{−1}), aimed at detecting Jupitermass companions (e.g., Bouchy et al. 2009; Moutou et al. 2014; Hébrard et al. 2016), as well as a search of smaller planets around Mdwarfs (e.g., Hobson et al. 2018, 2019; Díaz et al. 2019).
Since 2011, SOPHIE has been used for a survey of bright solartype stars, with the aim of detecting Neptunes and superEarths (Bouchy et al. 2011). For all the observations performed in this survey, the instrumental drift was measured and corrected for by recording on the detector, close to the stellar spectrum, the spectrum of a reference lamp. This one is a thoriumargon lamp before barycentric Jullian date (BJD) 2458181 and a FabryPerot interferometer after this date. The observations of HD 158259 were part of this program. Over the course of seven years, 290 measurements were obtained with an average error of 1.2 m s^{−1}. The data, which were corrected from instrumental drift and outliers (see Appendix A.1 and A.2), are shown in Fig. 1.
Fig. 1. SOPHIE radial velocity measurements of HD 158259 after outliers at BJD 2457941.5059, 2457944.4063, and 2457945.4585 have been removed. 
The RV is not the only data product that was extracted from the SOPHIE spectra. The SOPHIE pipeline also retrieves the bisector span (Queloz et al. 2001) as well as the (Noyes et al. 1984).
2.3. APT Photometry
Photometry has been obtained with the T11 telescope at the Automatic Photoelectric Telescopes (APTs), located at Fairborn Observatory in southern Arizona. The data, covering four observation seasons, are presented in more detail in Appendix A.3.
2.4. TESS results
HD 158259 was observed from 7 Oct to 2 Nov 2019 (sector 17). The TESS reduction pipeline (Jenkins et al. 2010, 2016) found evidence for a 2.1782 ±0.0006 d signal with a time of conjunction T_{c} = 2458766.049072 ± 003708 BJD (TOI 1462.01). The detection was made with a signaltonoise ratio of 8.05, which is above the detection threshold of 7.3 adopted in Sullivan et al. (2015).
3. Analysis of the data sets
3.1. Groundbased photometry and ancillary indicators
If the bisector span, , or photometry show signs of temporal correlation, in particular, periodic signatures, this might mean that the RVs are corrupted by stellar or instrumental effects (e.g., Queloz et al. 2001). To search for periodicities, we applied the generalized LombScargle periodogram iteratively (FerrazMello 1981; Zechmeister & Kürster 2009) as well as the ℓ_{1} periodogram (Hara et al. 2017) for comparison purposes. The process is presented in detail in Appendix A.4.
The periodogram presents a peak at 2900 d with a false alarm probability (FAP) of 4 × 10^{−12}. This signal is the only clear feature of the ancillary indicators, and it supports the presence of a magnetic cycle with a period ⩾1500 d. We note that both the APT photometry and bisector span present a peak around 11.6 d, though with a high FAP level. This periodicity as well as other periodicities found in the indicators were used to build candidate noise models for the analysis of the RVs, which is the object of the next section.
3.2. Radial velocities analysis
The RV time series we analyze here was corrected for instrument drift and from outliers. The process is described in Appendix A.1 and A.2. To search for potential periodicities, we computed the ℓ_{1} periodogram of the RV, as defined in Hara et al. (2017). This tool is based on a sparse recovery technique called the basis pursuit algorithm (Chen et al. 1998). The ℓ_{1} periodogram takes in a frequency grid and an assumed covariance matrix of the noise as input. It aims to find a representation of the RV time series as a sum of a small number of sinusoids whose frequencies are in the input grid. It outputs a figure which has a similar aspect as a regular periodogram, but with fewer peaks due to aliasing. The peaks can be assigned a FAP, whose intrepretation is close to the FAP of a regular periodogram peak.
The signals found to be statistically significant might vary from one noise model to another. To explore this aspect, we considered several candidate noise models based on the periodicities found in the ancillary indicators. The noise models were ranked with crossvalidation and Bayesian evidence approximations (see Appendix B). In Fig. 2, we show the ℓ_{1} periodogram of the SOPHIE RVs corresponding to the noise with the best cross validation score on a grid of equispaced frequencies between 0 and 0.7 cycle d^{−1}. The FAPs of the peaks pointed by red markers in Fig. 2 are given in Table 2. They suggest the presence of signals, in decreasing strengths of detection, at 366, 3.43, 7.95, 12.0, 5.19, 1920, 640 , and 17.4 days. The significance of the signals at 1.84, 17.7, and 34.5 d is found to be marginal to null.
Fig. 2. ℓ_{1} periodogram of the SOPHIE radial velocities of HD 158259 corrected from outliers (in blue). The periods at which the main peaks occur are represented in red. 
Periods appearing in the ℓ_{1} periodogram and their false alarm probabilities with their origin, the semiamplitude of corresponding signals, and M sin i with 68.7% intervals.
The FAPs reported in Table 2 were computed with a certain noise model, which is the best in a sense described in Appendix B. In this appendix, we also explore the sensitivity of the detections to the noise model choice. We find the detection of signals at 3.43, 5.19, 7.95, and 12.0 d to be robust. The detection of signals at 1920, 366, and 17.4 d is slightly less strong but still favored. There is evidence for a 640 d signal and hints of signals at 34.5 d and 1.84 d or 2.17 d. The latter two are aliases of one another. Indeed, the RV spectral window has a strong peak at the sidereal day (0.9972 d), which is common in RV time series (Dawson & Fabrycky 2010), and 1/1.840 + 1/2.177 = 1/0.9972 d^{−1}.
In Appendix C, the results are compared to a classical periodogram approach with a white noise model, which gives similar results but fails to unveil the 1.84/2.17 d candidate in the signal. We also study whether the apparent signals could originate from aliases of the periods considered here. This possibility is found to be unlikely.
We fit a model with a free error jitter and nine sinusoidal functions initialized at the periods listed in the upper part of Table 2. The data, which were phasefolded at the fitted periods, are shown in Figs. 3 and 4. The error bars correspond to the addition in quadrature of the nominal uncertainties and the fitted jitter.
Fig. 3. Radial velocity phasefolded at the periods of the signals appearing in the period analysis. From top to bottom: 2.17, 3.43, 5.19, 7.95, 12.0, and 17.4 d. 
Fig. 4. Radial velocity phasefolded, from top to bottom, at: 360, 750, and 2000 d. 
The residuals of the ninesignals model plus white noise fit have a root mean square (RMS) of 3.1 m s^{−1}, which is higher than the nominal uncertainties of the SOPHIE data (1.2 m s^{−1}). We studied the residuals with the methods of Hara et al. (2019) and found that they are temporally correlated, which, if not accounted for, might corrupt the orbital element estimates. In Appendix E, we show that a consistent model of the data can be obtained with a noise model that includes white and correlated components.
3.3. Periodicity origin
The origin of the 366, 640, and 1920 d signals is uncertain, and we do not claim planet detections at those periods. The 366 d signal is fully compatible with a yearly signal. Instrument systematics, such as the stitching effect (Dumusque et al. 2015), could produce this signal, and they are deemed to be its most likely origin. We fit a Gaussian process on the and used it as a linear predictor on the RVs, similarly to Haywood et al. (2014). When computing the periodogram on the residuals of the fit, there is no trace of signals in the 1500–3000 d and 600–800 d regions, so they might stem from a magnetic cycle. The signal at 34.5 d could be a faint trace of a planet near the 2:1 resonance, but its significance is too low to be conclusive.
The periods at 3.43, 5.19, 7.95, 12.0, and 17.4 d, which are significant in the RV analysis, most likely stem from planets. Here, we list five arguments that support this claim. (i) None of these periods clearly appear in the bisector span,, or photometry. (ii) While eccentric planets can be mistaken for planet pairs near a 2:1 MMR, they are very unlikely to appear as planets near a 3:2 resonance (Hara et al. 2019). (iii) The periods could be due to instrument systematics. We find it unlikely since the periods of the planet candidates do not consistently appear in the 123 other data sets of the survey HD 158259 is a part of Hara et al. (in prep.). (iv) All the signals are consistent in phase and amplitude (see Appendix D). (v) Most importantly, the period ratio of two subsequent planets is very close to 3:2, namely 1.51, 1.53, 1.51, and 1.44, and they have very similar estimated masses (see Sect. 3.4). Pairs of planets close to the 3:2 period ratio are known to be common (Lissauer et al. 2011a; Steffen & Hwang 2015), and it seems unlikely that the stellar features would mimic this specific spacing of periods.
Signals at 12 and 17.4 d verify the five points listed above and are thus considered as planets. We, however, point out that the predicted rotation period of the star is 18 ± 5 d. This means that signals could be present in this period range, not necessarily at the stellar rotation period or its harmonic (Nava et al. 2020). Furthermore, in Appendix B, we show that the 17.4 d signal is less significant and that certain noise models favor 17.7 d over 17.4 d. It appears that fitting signals at 17.4 or 17.7 d does not completely remove the other. This might point to differential rotation or dynamical effects, but it is most likely due to modeling uncertainties. Nonetheless, the nature and period of the 17.4 d signal are subject to a little more caution than the other planets. It is thus conservatively classified as a strong planet candidate.
The 2.17 d signal appearing in TESS has a counterpart in the RVs, which appear at 1.84 d (alias) here. The RV signal is marginally significant but compatible with the observation of a transit. Indeed, the time of conjunction as measured by TESS is BJD 2458766.049072 ± 0.003708. At this epoch, the mean longitude of the innermost planet predicted from the RV with a prior on the period set from TESS data is (see Sect. 3.4 for details on the posterior calculation). On this basis, the 2.17 d signal is considered to stem from a planet. We note that no trace of the other planets is found in the TESS data.
In Fig. 5 we represent the configuration of the system at the time of conjunction given by the TESS pipeline. The markers represent the position of the planets corresponding to their posterior mean (semi major axis and mean longitude). The marker size is proportional to an estimated radius ∝(m sin i)^{0.55}, following Bashi et al. (2017). The plain colored lines correspond to 68% credible intervals on the mean longitude.
Fig. 5. Configuration of the system estimated from the RV at the time of conjunction estimated by TESS (BJD 58766.049072). 
In summary, we deem six of the nine significant signals to originate from planets: 2.17, 3.43, 5.19, 7.95, 12.0, and 17.4 d, with a slightly lower confidence in 17.4 d. From now on, they are referred to as planets b, c, d, e, f, and (strong) candidate g, respectively.
3.4. Orbital elements
To derive the uncertainties on the orbital elements of the planets, we computed their posterior distribution with a Monte Carlo Markov chain algorithm (MCMC). The model includes signals at 2.17, 3.43, 5.19, 7.95, 12.0, and 17.4 d, a linear predictor fitted on the with a Gaussian process analysis, as well as a correlated noise model with an exponential decay. The prior on eccentricity was chosen to strongly disfavor e > 0.1 (see Sect. 4 for justification). The details of the posterior calculations are presented in Appendix E and the main features of the signals are reported in Table 2. Planets c, d, e, f, and planet candidate g exhibit m sin i ≈ 6 M_{⊕} and planet b has m ≈ 2 M_{⊕} and R ≈ 1.2 R_{⊕}, which is compatible with a terrestrial density. The yearly signal has an amplitude of ≈3 m s^{−1}. We find a significantly nonzero timescale of the noise of d. A more detailed analysis shows hints of nonstationary noises (see Appendix D).
The TESS photometry exhibits a transit signal from b, but not from c, d, e, f, or g. Given the masses of c, d, e, f, and g, their transits should have been detected had they occurred. Assuming a coplanar system, this can be explained if the inclination departs from 90° from at least (so that the transit of c cannot be detected) to at most (the transit of b must be seen). This would translate to a sin i between 0.985 and 0.993. Alternately, the planets could be relatively inclined.
4. Dynamical analysis
4.1. Stability
For the dynamical analysis of the system, we consider the planets b, c, d, e, f, and candidate g, whose period ratios are close to the 3:2 MMRs: P_{c}/P_{b} = 1.57, P_{d}/P_{c} = 1.51, P_{e}/P_{d} = 1.53, P_{f}/P_{e} = 1.51, and P_{g}/P_{f} = 1.44. To check the existence of stable solutions, we proceeded as follows. We performed the MCMC analysis with exactly the same priors as in Sect. 3.4 except with a looser prior on eccentricity. Each MCMC sample was taken as an initial condition for the system. Its evolution was integrated 1 kyr in the future using the 15th order Nbody integrator IAS15 (Rein & Spiegel 2015) from the package REBOUND^{1} (Rein & Liu 2012). General relativity was included via REBOUNDx, using the model of Anderson et al. (1975). In total, 36 782 system configurations were integrated. We consider the system to be unstable if any two planets have an encounter below their mutual Hill radius. A total of 1700 samples lead to integrations where the stability condition was not violated. We computed their empirical probability distribution functions and found constraints on the eccentricities ≲0.1.
4.2. Resonances
The planets of HD 158259 could be locked in 3:2 or three body resonances. This would translate to the following conditions. We consider three subsequent planets, indexed by i = 1, 2, 3 from innermost to outermost. With λ_{i}, P_{i}, and ϖ_{i}, we denote their mean longitudes, periods, and arguments of periastron. In case of two body resonances, the angle
would librate. Following Delisle (2017), in case of threebody resonances, the socalled Laplace angle
should librate around π. The posterior distributions of ψ_{ij} and ϕ_{ijk} as well as their derivatives were computed for any doublet and triplet of planets from the MCMC samples, assuming that ϖ_{j} is constant on the timescale of the observations. We conclude that the system is not locked in two and three body resonances. We note that the circulation of the angles occurs at different rates, ϕ_{cde} being the slowest.
Nevertheless, period ratios so close to 3:2 are very unlikely to stem from pure randomness. It is therefore probable that the planets underwent migration in the protoplanetary disk, during which each consecutive pair of planets was locked in 3:2 MMR. The observed departure of the ratio of periods of two subsequent planets from exact commensurability might be explained by tidal dissipation, as was already proposed for similar Kepler systems (e.g., Delisle & Laskar 2014). Stellar and planet mass changes have also been suggested as a possible cause of resonance breaking (Matsumoto & Ogihara 2020). The reasons behind the absence of threebody resonances, which are seen in other resembling systems (e.g., Kepler80, MacDonald et al. 2016), are to be explored.
5. Conclusion
In this work, we have analyzed 290 SOPHIE measurements of HD 158259. The analysis of the radial velocity data, including a correction of the instrument drift and over 7500 correlated noise models, support the detection of four planets (c, d, e, and f) and a strong candidate (g), with respective periods of 3.43, 5.19, 7.95, 12.0, and 17.4 d. They all exhibit a ≈6 M_{⊕} mass. There is substantial evidence for a planetary origin of the 17.4 d signal, and the remaining concerns on its nature should be cleared with a better estimate of the stellar period. Furthermore, the TESS data exhibit a 2.17 d signal that is compatible with the RVs, which allows one to claim the detection of an additional planet (b) and to measure a density of ρ_{⊕}. The TESS photometry does not show other transits. There exist stable configurations of the sixplanet system that are compatible with the error bars.
While many compact nearresonance chains have been detected by transits, they have been rare in radial velocity surveys so far. The present analysis shows that they can be detected, provided there are enough data points and an appropriate accounting of correlated noises (instrumental and stellar).
HD 158259 b, c, d, e, f, and g are such that subsequent planets have period ratios of 1.57, 1.51, 1.53, 1.51, and 1.44 with increasing period. Subsequent planet pairs and triplets are close to, but not within, 3:2 and threebody mean motion resonances. The period ratios are consistent with the distribution of period ratios of planet pairs found in Kepler, which exhibits a peak at 1.52 (see Lissauer et al. 2011a; Fabrycky et al. 2014; Steffen & Hwang 2015). The similarity of the masses of the planets of the system is compatible with the hypothesis that planets within the same system have similar sizes (Lissauer et al. 2011b; Ciardi et al. 2013; Millholland et al. 2017; Weiss et al. 2018). The configuration of the planets can be explained by existing formation scenarios (e.g., Terquem & Papaloizou 2007; Delisle & Laskar 2014). The proximity to 3:2 resonances of the HD 158259 system is reminiscent of Kepler80 (MacDonald et al. 2016). However, the latter presents threebody resonances, while HD 158259 does not. It would be interesting to investigate scenarios that can explain the differences between the two systems.
The REBOUND code is freely available at http://github.com/hannorein/rebound.
Acknowledgments
We warmly thank the OHP staff for their support on the observations. N.C.H. and J.B. D. acknowledge the financial support of the National Centre for Competence in Research PlanetS of the Swiss National Science Foundation (SNSF). This work was supported by the Programme National de Planétologie (PNP) of CNRS/INSU, cofunded by CNE. G.W.H. acknowledges longterm support from NASA, NSF, Tennessee State University, and the State of Tennessee through its Centers of Excellence program. X.De., X.B., I.B., and T.F. received funding from the French Programme National de Physique Stellaire (PNPS) and the Programme National de Planétologie (PNP) of CNRS (INSU). X.B. acknowledges funding from the European Research Council under the ERC Grant Agreement n. 337591ExTrA. This work has been supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). This work is also supported by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR15IDEX02), through the funding of the ‘Origin of Life’ project of the Univ. GrenobleAlpes. V.B. acknowledges support from the Swiss National Science Foundation (SNSF) in the frame of the National Centre for Competence in Research PlanetS, and has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (project Four Aces; grant agreement No 724427). This work was supported by FCT – Fundação para a Ciência e a Tecnologia/MCTES through national funds (PIDDAC) by this grant UID/FIS/04434/2019, as well as through national funds (PTDC/FISAST/28953/2017 and PTDC/FISAST/32113/2017) and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalizaȩão (POCI010145FEDER028953 and POCI010145FEDER032113). N.AD. acknowledges support from FONDECYT #3180063. X.Du. is grateful to the Branco Weiss Fellowship–Society in Science for continuous support.
References
 Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Amelunxen, D., Lotz, M., McCoy, M. B., & Tropp, J. A. 2014, Inf. Inference J. IMA, 3, 224 [CrossRef] [Google Scholar]
 Anderson, J. D., Esposito, P. B., Martin, W., Thornton, C. L., & Muhleman, D. O. 1975, ApJ, 200, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Baluev, R. V. 2008, MNRAS, 385, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Bashi, D., Helled, R., Zucker, S., & Mordasini, C. 2017, A&A, 604, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Eggenberger, A., Santos, N. C., et al. 2010, A&A, 523, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Bouchy, F., Hébrard, G., Delfosse, X., et al. 2011, EPSCDPS Joint Meeting 2011, 2011, 240 [NASA ADS] [Google Scholar]
 Bouchy, F., Hebrard, G., Udry, S., Delfosse, X., & Boisse, I. 2009, A&A, 505, 853 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cannon, A. J., & Pickering, E. C. 1993, VizieR Online Data Catalog: III/135A [Google Scholar]
 Chandler, C. O., McDonald, I., & Kane, S. R. 2016, VizieR Online Data Catalog: J/AJ/151/59 [Google Scholar]
 Chen, S. S., Donoho, D. L., & Saunders, M. A. 1998, SIAM J. Sci. Comput., 20, 33 [Google Scholar]
 Ciardi, D. R., Fabrycky, D. C., Ford, E. B., et al. 2013, ApJ, 763, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Courcol, B., Bouchy, F., Pepe, F., et al. 2015, A&A, 581, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Delisle, J.B. 2017, A&A, 605, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., & Laskar, J. 2014, A&A, 570, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Ségransan, D., Dumusque, X., et al. 2018, A&A, 614, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J. B., Hara, N. C., & Ségransan, D. 2019, A&A, submitted [Google Scholar]
 Delisle, J. B., Hara, N., & Ségransan, D. 2020, A&A, 635, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Díaz, R. F., Delfosse, X., Hobson, M. J., et al. 2019, A&A, 625, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., & Latham, D. W. 2015, ApJ, 808, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [NASA ADS] [CrossRef] [Google Scholar]
 FerrazMello, S. 1981, AJ, 86, 619 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration 2018, VizieR Online Data Catalog: I/345 [Google Scholar]
 Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., & Correia, A. C. M. 2017, MNRAS, 464, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738 [NASA ADS] [CrossRef] [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [NASA ADS] [CrossRef] [Google Scholar]
 Hébrard, G., Arnold, L., Forveille, T., et al. 2016, A&A, 588, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henry, G. W. 1999, PASP, 111, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Hobson, M. 2019, Exoplanet detection around M dwarfs with near infrared and visible spectroscopy [Google Scholar]
 Hobson, M. J., Díaz, R. F., Delfosse, X., et al. 2018, A&A, 618, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hobson, M. J., Delfosse, X., AstudilloDefru, N., et al. 2019, A&A, 625, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 Jenkins, J. M., Chandrasekaran, H., McCauliff, S. D., et al. 2010, in Proc. SPIE, 7740, 77400D [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Software and Cyberinfrastructure for Astronomy IV, Proc. SPIE. 9913, 99133E [Google Scholar]
 Jones, D. E., Stenning, D. C., Ford, E. B., et al. 2017, ArXiv eprints [arXiv:1711.01318] [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [Google Scholar]
 Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011a, ApJS, 197, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011b, ApJS, 197, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, T. A., Barros, S. C. C., Santerne, A., et al. 2019, A&A, 631, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
 MacDonald, M. G., Ragozzine, D., Fabrycky, D. C., et al. 2016, AJ, 152, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, Y., & Ogihara, M. 2020, ApJ, accepted [arXiv:2003.01965] [Google Scholar]
 Mayor, M., Udry, S., Lovis, C., et al. 2009, A&A, 493, 639 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Mortier, A., & Collier Cameron, A. 2017, A&A, 601, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moutou, C., Hébrard, G., Bouchy, F., et al. 2014, A&A, 563, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nava, C., LópezMorales, M., Haywood, R. D., & Giles, H. A. C. 2020, AJ, 159, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, B. E., Ford, E. B., Buchner, J., et al. 2020, AJ, 159, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Noyes, R. W. 1984, in Space Research in Stellar Activity and Variability, eds. A. Mangeney, & F. Praderie, 113 [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning Adaptive Computation and Machine Learning (The MIT Press) [Google Scholar]
 Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, J. Astron. Telescopes Instrum. Syst., 1, 1 [NASA ADS] [Google Scholar]
 Schuster, A. 1898, Terr. Mag., 3, 13 [Google Scholar]
 Schwarz, G. 1978, Ann. Statist., 6, 461 [Google Scholar]
 Shallue, C. J., & Vanderburg, A. 2018, AJ, 155, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, S. S., & Wilk, M. B. 1965, Biometrika, 52, 591 [Google Scholar]
 Steffen, J. H., & Hwang, J. A. 2015, MNRAS, 448, 1956 [NASA ADS] [CrossRef] [Google Scholar]
 Sullivan, P. W., Winn, J. N., BertaThompson, Z. K., et al. 2015, ApJ, 809, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 [NASA ADS] [CrossRef] [Google Scholar]
 Tibshirani, R. 1994, J. R. Stat. Soc. Ser. B, 58, 267 [Google Scholar]
 Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, J.W. 2013, ApJS, 208, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Complementary analysis of the time series
A.1. Accounting for instrumental effects in RV
SOPHIE experiences a drift of the zero velocity point due to several factors, as follows: a change of fiber, calibration lamp aging, and other systematic effects. A drift estimate was obtained by observing reference stars, which are deemed to have a nearly constant velocity, each night of observations. The reference stars’ velocities were combined and interpolated to create an estimate of the drift as a function of time. This one was then subtracted from all the time series of the observation program. The estimation procedure, which is similar to that of Courcol et al. (2015), is presented in detail in Hara et al. (in prep.). We performed the timeseries analysis of the data obtained with the original correction of Courcol et al. (2015). The results are very similar; the only notable difference is a lower statistical significance of the 12.0 days signal.
A.2. Data selection
Some of the data points of the radial velocity and time series are excluded from the analysis. In this section, we present the methodology adopted.
We removed outliers from the time series with a criterion based on the median absolute deviation (MAD). We computed σ = 1.48 MAD, which is the relation between the standard deviation σ and the MAD of a Gaussian distribution. We excluded the data points if their absolute difference to the median is greater than k × σ with k = 4.
In Fig. A.1, we show the radial velocity data points and their error bars based on the nominal uncertainties of the measurements. The three points excluded from the analysis based on the MAD criterion are represented in green. There is one point with a very clear deviation and two outliers with a lower deviation. The farthest outlier prevents finding any significant signal in the RVs besides a signal with a timescale of 2000 d. Including the two other outliers in the RV analysis yields minor changes in the analysis, and it does not challenge the detection of the planets. The and bisector measurements at the dates of the three outliers are also excluded from the analysis.
Fig. A.1. SOPHIE radial velocities and nominal 1σ error bars, the three points which were removed are shown in green. 
In the case of the , besides the removal criterion based on the MAD, the measured after BJD 2458000 are excluded. After this date, the values of the are not reliable. This is due to the change of the calibration lamp from thoriumargon to FabryPerot, which leaks on the stellar spectrum as a continuum background. No effects are seen in the RV of the constant stars, which were observed each night, but the leaking affects the measurements of the . This issue is described in greater detail in Hobson (2019), p. 92. We note that excluding some of the points has a very limited impact on our results, which do not depend on the exact modeling of low frequency components in the RV signal. Furthermore, excluding the RV points after the calibration change does not affect our results. The points selected for the analysis are presented in Fig. A.2 in blue, while the points excluded are represented in red. In our analysis, we do not exclude the RV points at the dates of the points excluded from the . The application of the MAD criterion on the bisector span and photometry does not lead to the removal of any point.
Fig. A.2. measurements of HD 158259. Points in blue are conserved for the analysis; the points in red are excluded from the analysis. 
A.3. APT Photometry
We acquired 320 observations of HD 158259 during the 2002, 2003, 2004, and 2005 observing seasons with the T11 0.80 m Automatic Photoelectric Telescope (APT) at Fairborn Observatory in southern Arizona. T11 is equipped with a twochannel precision photometer that uses a dichroic mirror and two EMI 9124QB bialkali photomultiplier tubes to measure the Strömgren b and y passbands simultaneously. We observed the program star with respect to three comparison stars and computed the differential magnitudes as the difference in brightness between the program star and the mean brightness of the two best comparison stars. To improve the photometric precision, we combined the differential b and y magnitudes into a single (b + y)/2 passband. The precision of a single observation is typically around 0.0015 mag. The T11 APT is functionally identical to our T8 0.80 m APT described in Henry (1999), where further details of the telescope, precision photometer, and observing and data reduction procedures can be found.
Table A.1 gives a summary of the photometric results. No significant variability is found within any observing season. The seasonal means exhibit a range of 0.0018 mag, which is probably due to similar variability seen in the mean magnitudes of the comparison stars. The lack of observed spot activity is consistent with the star’s low value of . For further analysis, the four observing seasons of APT photometry were normalized such that the last three seasons have the same mean magnitude as the first (see next section).
Summary of APT photometric observations for HD 15825.
A.4. Period search
The presence of signals in the ancillary indicators might give hints as to the instrumental and stellar features in the radial velocity. We are particularly interested in periodic signals, which could also appear in the RV and mimic a planetary signal. To search for periodic signals, we analyzed the photometry and ancillary indicators with classical and ℓ_{1} periodograms. The periods found are used to build candidate quasiperiodic noise models in the RV.
We first present the classical periodogram analysis. We computed the generalized LombScargle periodogram (FerrazMello 1981; Zechmeister & Kürster 2009) of the ancillary indicators on a grid of frequencies spanning from 0 to 1.5 cycles day^{−1}. We report the strongest periodic signatures. The false alarm probability (FAP) of the highest peaks of the periodograms were computed using the Baluev (2008) analytical formula. After computing the periodogram and checking that no significant high frequency signal is found, we performed the search on a grid of frequencies from 0 to 0.95 cycles day^{−1} to avoid aliases in the one day region. We subtracted signals at the periods found iteratively and computed the periodograms of the residuals. We also applied ℓ_{1} periodograms for comparison purposes.
The analysis was performed after the data selection, which is presented in Appendix A.2. The periodogram (see Fig. A.5) presents a clear longterm trend and, besides peaks in the one day region, it peaks at 119 and 64 d. The iterative period search on a frequency grid from 0 to 0.95 cycle day^{−1} gives signals at 3500, 119, and 32 d, with FAPs 4 × 10^{−12}, 0.18 and 0.14. The longterm signal in the might correspond to a magnetic cycle. In Sect. B.2, we fit a Gaussian process model to the data (see Fig. B.6), which was used in the RV processing.
The bisector span periodogram is presented in Fig. A.4. It presents peaks, in order of decreasing amplitude at 0.9857, 552, 85, 1576, and 23.5 days, where 0.9857 is an alias of 85 days. The FAP of the highest peak is 0.08, which indicates marginal evidence against the hypothesis that the bisector behaves like white noise. The iterative period search from 0 to 0.95 cycle day^{−1} points to signals at 550, 85, and 11.5 d with FAPs of 0.13, 0.21, and 0.25.
In Fig. A.3, we represent the periodogram of the photometric data. The maximum peak occurs at a period of 0.979 days, which is an alias of 11.6 days, and has a FAP of 0.33. The iterative search yields 11.63 and 108 days as dominant periods with FAPs of 0.75 and 0.43. The fact that 11.5 days both appear in the bisector and photometry might indicate that there is a weak stellar feature at this period, which is possibly due to the rotation period or one of its harmonics. The fitted amplitudes of sine functions initialized at the orbital periods are given in Table A.2. In each case, the peaktopeak amplitude is very small (a fraction of a millimagnitude) and of the same order as its uncertainty.
Fig. A.3. Periodogram of the photometric data computed between on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 0.979 d (in red). The four subsequent tallest peaks are, in decreasing order, at 0.817 11.636, 1.028, and 108 d (in yellow). 
Fig. A.4. Periodogram of the bisector span computed on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 0.985 d (in red). The four subsequent tallest peaks are, in decreasing order, at 552, 85.5, 1576, and 23.48 d (in yellow). 
Fig. A.5. Periodogram of the computed on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 2900 d (in red). Besides aliases in the one day region, there is a peak at 119 and 64 days. 
Photometric amplitudes on the planetary orbital periods.
The ℓ_{1} periodograms of the photometry, bisector span, and the (resp. Figs. A.6–A.8) were computed on a frequency grid from 0 to 0.95 cycles day^{−1} to avoid the one day region. The ℓ_{1} periodogram aims to express the signal as a sum of a small number of sinusoidal components. Here, the ℓ_{1} periodograms present numerous peaks, which is characteristic of noisy signals, as they do not have a sparse representation in frequency. The periods appearing in ℓ_{1} periodograms are consistent with those appearing in the classical analysis. We simply note the addition of a signal at 23.5 d in the bisector span ≈11.6 × 2.
Fig. A.6. ℓ_{1} periodogram of the APT photometric measurements. 
Fig. A.7. ℓ_{1} periodogram of the bisector span. 
Fig. A.8. ℓ_{1} periodogram of the . 
In conclusion, the analysis of the photometry and ancillary indicators supports the existence of a longterm magnetic cycle, appearing in the periodogram, with a period ≥1500 d. The period cannot be resolved due to the timespan of the SOPHIE observations (2560 days). The presence of several marginally significant periods in the bisector span might indicate the presence of correlated noise in the bisector, and possibly in the RV.
Appendix B: Impact of the noise model on the planet detection
B.1. Selection with cross validation
In this section, we present the noise model chosen and the sensitivity of the detection to the noise model. The various sources of noise in the RV were modeled as in Haywood et al. (2014) by a correlated Gaussian noise model. In practice, one chooses a parametrization of the covariance matrix of the noise V(θ), where the element of V at index k, l depends on t_{k} − t_{l} and a vector of parameters θ. In the following analysis, the parametrization chosen for V is such that its element at index k, l is
where σ_{k} is the nominal measurement uncertainty, σ_{W} is an additional white noise, σ_{C} is a calibration noise, and where c equals one if measurements k and l are taken within the same night and zero otherwise. The quantities σ_{R} and τ_{R} parametrize a correlated noise, which might originate from the star or the instrument. Additionally, σ_{QP}, τ_{QP}, and P_{act} parametrize a quasiperiodic covariance, potentially resulting from spots or faculae. The form of this covariance is compatible with the spleaf software (Delisle et al. 2020) and, except for the calibration noise, the CELERITE model (ForemanMackey et al. 2017).
To study the sensitivity of the detection to the noise model θ = (σ_{W}, σ_{C}, σ_{R} τ_{R}, σ_{QP}, τ_{QP}, P_{act}), we proceeded as follows. We first considered a grid of possible values for each component of θ. For instance, σ_{R} = 0, 1, 2, 3, 4 m s^{−1}, τ_{R} = 0, 1, 6 d etc. The θ with all the possible combinations of the values of its components were generated, and the corresponding covariance matrices were created according to Eq. (B.1).
The particular values taken in the present analysis are reported in Table B.1. The decay time scales of the red and quasi periodic components (subscripts R and QP) include 0 in which case they are white noise jitters. The σ_{R} and σ_{QP} were chosen such that when τ_{R} = τ_{QP} = 0, there exists a value of that is greater than the total variance of the data, and we subdivided the possible values of σ_{R}, σ_{QP} in smaller steps. The correlation timescales of τ_{R} = 0, 1, and 6 d correspond to no correlation, daily correlation, and noise correlated on a whole run of observations, which is typically 6 days. The P_{act} candidate corresponds to peaks that appear in the period analysis of the radial velocity or ancillary indicators (Sect. A.4). We note that 2000 d was not considered as it is degenerate with an exponential decay due to the observation time span.
Value of the parameters used to define the grid of models tested.
For each matrix in the list of candidate noise models, the ℓ_{1} periodogram was computed, and the frequencies that have a peak with FAP < 0.05 were selected. We then attributed a score to the couple (noise model, planets) with two different metrics.
The first metric is based on cross validation. We selected 70% of the data points randomly – the training set – and fit a sinusoidal model at the selected frequencies on this point. On the remaining 30% – the test set – we computed the likelihood of the data knowing the model fitted. The operation of selecting of training set randomly and evaluating the likelihood on the test set was repeated 200 times. We took the median of the 200 values of the likelihood as the cross validation score of the noise model. As a result of this procedure, each noise model has a cross validation score (CV score).
The second metric is the approximation of the Bayesian evidence of the (noise model, planets) couple. For a given covariance matrix, we fit a sinusoidal model initialized at the periods selected by the FAP criterion. We then made a second order Taylor expansion of the posterior (with priors set to one in all the variables) and estimated the evidence of the model, which is also called the Laplace approximation (Kass & Raftery 1995). The procedure is explained in greater detail in Nelson et al. (2020), Appendix A.4.
In Fig. B.2, we represent the histogram of the values of the CV score for all the noise models considered. The models with the 20% highest CV score are represented in blue, and they present similar values of the CV score. We call the set of these models CV_{20}.
Each model of CV_{20} might lead to different peaks selected with the FAP < 0.05 criterion. The periods and FAPs of the peaks selected are represented in Fig. B.3 by the yellow points. For a comparison with the best models, that is, the model with the maximal CV score, the periods appearing in the best model (red dots in Fig. 2) are represented by the blue dashed lines. We represent the median of the FAPs for each of these periods in the CV_{20} and their FAP in the best models (purple and red, respectively). These values are also listed in Table B.2, where we also give the percentage of cases where a period corresponds to a FAP < 0.05 in CV_{20}.
Periods appearing in the ℓ_{1} periodogram of the model with the highest CV score and their false alarm probabilities.
We plotted the same Fig. B.3 and the same Table B.2 for the 20% best models in terms of evidence (E_{20}) (Fig. B.4 and Table B.3). The ℓ_{1} periodogram corresponding to the highest ranked model is shown in Fig. B.1). Peaks at 3.4, 5.2, 7.9, and 12 d have similar behaviors with cross validation and approximate evidence. We see three notable differences when ranking a model with evidence: ≈10% of the selected models display significant peaks at 1.84 or 2.17 d, 17.4 reaches a 5% FAP in only 38% of E_{20}, and 8% of the models favor 17.4 over 17.7 d. Finally, the inclusion percentage and average FAP of long period signals significantly drops. We have also ranked models with AIC (Akaike 1974) and BIC (Schwarz 1978), which yield results very similar to approximate evidence.
Fig. B.1. ℓ_{1} periodogram corresponding to the highest approximate evidence. 
Fig. B.2. Histogram of the values of the cross validation score. Best 20% and lowest 80% are represented in blue and orange, respectively. 
Fig. B.3. FAPs of the peaks of the 20% best models in terms of cross validation score that have a FAP > 0.05. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 
Fig. B.4. FAPs of the peaks of the 20% best models in terms of Laplace approximation of the evidence that have a FAP > 0.05. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 
Periods appearing in the ℓ_{1} periodogram of the model with the highest approximated evidence and their false alarm probabilities.
When using a LASSOtype estimator (Tibshirani 1994), such as the ℓ_{1} periodogram, for a fixed dictionary (here, a fixed noise model), it is common practice to select the model with crossvalidation as the solution follows the socalled LASSO path. Here, this would constitute a viable alternative to selecting the peaks that have FAP < 0.05. This was tested on a grid of parameters such as B.1, and it yields very similar conclusions.
B.2. Longterm model
It has been found that the has a strong longterm signal. This one can be included in the analysis with a Gaussian process analysis, similarly to Haywood et al. (2014). We consider a simple covariance model for the
We assume that σ_{R} is equal to the standard deviation of the data, and we fit the parameters σ_{W} and τ. We find σ_{W} = 0.9σ_{R} and τ = 770 d. We then predicted the Gaussian process value and its covariance with formulae 2.23 and 2.24 in Rasmussen & Williams (2005). In Fig. B.6, we represent the raw data on which the fitting is made. The Gaussian process and the one sigma standard deviation of the marginal distribution at t are represented in light blue, and the prediction at the radial velocity measurement times is in orange.
We then performed the same analysis as in Sect. B.1, except that we include the smoothed (orange points in Fig. B.6) as a linear predictor in the model. The results are presented in Fig. B.5. These are almost identical to the results of Sect. B.1, except that the 2000 d and 640 d signals disappear.
Fig. B.5. FAPs of the peaks of the 20% best models that have a FAP > 0.05, when adding a linear activity model fitted as a Gaussian process. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 
Fig. B.6. Raw and smoothed time series. The dark blue points represent the raw data used for the prediction. The light blue lines represent the Gaussian process prediction and its ±1σ error bars (see Eq. (B.2)). The orange points are the predicted values of the Gaussian process at the radial velocity measurement times. 
B.3. Discussion
The noise model has a strong impact on the false alarm probabilities. When instantiating the covariance matrix with an exponential kernel, such as Eq. (B.2), the characteristic time τ can be seen as setting the threshold between “short” and “long” periods. The FAPs of the signals with periods lower than τ are lower and higher, respectively, compared to the ones computed with a white noise model. Using a correlated noise model reduces the chances of spurious detection claims at long periods, but on the other hand, it decreases the statistical power (see Delisle et al. 2020 for a more detailed discussion on that point). The procedure of Sect. B.1 aims to balance the two types of errors: If the white noise model cannot be excluded (here, this means being in the best 20% noise models), then the short and long period signal appear with lower and higher FAPs, respectively, resulting in a compromise value. We have found the result to be insensitive to the exact threshold on the models selected (choosing the best 10, 20, and 30% yields identical conclusions). The interpretation of the differences between cross validation and evidence is that evidence appears to favor models with stronger correlated components, which damp the amplitude of long period signals and boost high frequency signals.
We have found that signals at 3.43, 5.19, 7.95, 12.0, 17.4, and 366 d consistently appear in the models. The 1920 and 640 d signals appearing in Sect. B.1 disappear when modeling the activity as in Sect. B.2, pointing to a stellar origin. The strength of the 17.4 d signal varies with the model ranking technique, which makes its detection less strong than the other signals. We, however, note that this one was put in competition with a quasiperiodic model at the same period, which is not the case of the other signals. Overall, we claim that the 3.43, 5.19, 7.95, 12.0, 17.4, and 366 d periodicities are present in the signal and there is a candidate at 1.84/2.17 d.
As a remark, the noise model selection procedure is close to Jones et al. (2017), which ranks noise models with cross validation, BIC, and AIC. The difference lies in the fact that instead of comparing noise models with free parameters, we compare models that are couples of the noise model with fixed parameters and planets with fixed periods.
Appendix C: Periodogram analysis
C.1. Iterative search
In this section, we perform the period search in an iterative way. We follow a standard planet search algorithm. We first compute at which frequencies the two maximum peaks of the spectral windows are attained (besides the peak in ω = 0) and denote them by ω_{S1} and ω_{S2}. Here, ω_{S1} = 1.0027 and ω_{S2} = 0.999 cycles day^{−1}.
The planet count is initialized at 0. We compute the generalized LombScargle periodogram, and we search at which frequency its maximum is attained ω_{0}. We then fit a Keplerian signal initialized at ω_{0}, but depending on the value of ω_{0}, we might choose to fit ω_{0} − ω_{Si} instead, where i = 1, 2. We then add one to the planet count and compute the periodogram of the residuals. The frequency of the maximum peak is denoted by ω_{1}; we then fit two Keplerian functions initialized at ω_{0} and ω_{1}, compute the periodogram of the residuals, and so on until a non significant signal is found.
In Fig. C.1 we represent the subsequent periodograms, on a period range spanning from 0 to 1.5 cycles day^{−1}, which were computed iteratively. The periods selected at each iteration are marked in red, and their false alarm probability corresponding to their peaks is given. It might happen that the maximum of the periodogram occurs at a period near 1 day, in which case the corresponding period is highlighted by a vertical yellow line. Since in each occurrence of this situation, there are peaks at a longer period which are aliases of these ones, the peaks at longer periods are preferred.
Fig. C.1. Subsequent periodograms after a fit of a Keplerian orbit model initialized at the maximum peak of the periodogram 
Overall, the periodogram results are similar to those of the ℓ_{1} periodogram, the most striking difference being that in the classical approach, the 7.95 dsignal does not appear, while the 17.4 dsignal has a stronger significance. Furthermore, the 1.84/2.17 candidate does not appear. In the classical approach, the 7.95 dsignal is absorbed in the Keplerian fit initialized at 17.4 d. When performing the iterative search with sinusoids only, the 7.95 dsignal appears at the fifth iteration. The decreased significance of the 17.4 dperiod in the ℓ_{1} periodogram compared to the classical periodogram stems from accounting for correlated noises in the ℓ_{1} periodogram. As discussed in Sect. B.3, the signals with a period greater than the timescale of the noise see their significance decrease with respect to a white noise model.
C.2. Aliases
In the previous section, we have seen that at the fourth iteration, the highest peak in the periodogram occurs at 1.238 d, which, given the spectral window, is an alias of 5.198 d. Furthermore, when performing an iterative search with a circular model between 0 and 1.5 cycles day^{−1}, the subsequent maximum periodogram peaks include aliases of 3.4, 7.95, 17.4, and 2000 d. This raises the question of whether the periodicities claimed here might in fact originate from shorter periods.
It seems improbable that any of the planets would be at the aliases close to one day. Indeed, planet couples with period ratios near 1.52 are found to be common (Lissauer et al. 2011a; Fabrycky et al. 2014; Steffen & Hwang 2015). If the chain found was spurious, the aliases would have to fall by chance such that a 3:2 chain would be compatible with the data. This cannot be completely excluded, but it seems improbable.
Appendix D: Phase and amplitude consistency
A property of planetary signals is that, aside from signatures of gravitational interaction between planets, their amplitude and phase do not vary with time. To check the consistency in phase and the amplitude of the signals found, we performed two tests.
First, we computed the ℓ_{1} periodogram of the signal starting with the first 60 points and adding 1 point at a time up to 287, with the noise model corresponding to the best cross validation score, similarly to the procedures suggested in Schuster (1898) and Mortier & Collier Cameron (2017). From ≈180 points, the signals appearing in the ℓ_{1} periodogram are identical to those from Table 2. This result is also consistent when it is performed backwards, starting with the last 60 points. In Fig. D.1, we represent the evolution of the peak amplitudes at several periods as a function of the number of points n in the vicinity of periods of interest. By vicinity, we mean that the value plotted is the maximum of the peak occurring within 1/T_{obs(n)} in frequency of the period of interest; T_{obs(n)} being the total observation time of the data set including n points. It appears that the amplitude of the signals reported as planets and planetary candidates increases steadily (top panel), while other signals show more variability (bottom panel). We note that before 180 points, the frequency resolution is insufficient in distinguishing signals at 640 and 1920 d.
Fig. D.1. Evolution of the ℓ_{1} periodogram amplitudes of peaks at different periods. Top: planets and planet candidates. Bottom: other signals. The black dashed line indicates the transition from calibration with thoriumargon to FabryPerot. 
It appears that most of the candidates appear from ≈130 points. This could be due to a phenomenon analogous to phase transitions observed in sparse recovery with random matrices, where the number of signals appropriately recovered changes sharply in the vicinity of a critical number of observations (Amelunxen et al. 2014). However, the ℓ_{1}periodogram of the first, middle, and last 150 points show important differences. The first 150 points fail to unveil most of the signals that are later confirmed, while they are seen in the ℓ_{1}periodogram of the middle and last 150 points.
We now check the consistency in phase of the signals. We consider the points up to BJD 2457529.4867 (excluded) and after this date, so that we have two data sets of 144 and 143 points. The model described in Sect. 3.4 was fitted onto each half of the data. The likelihood and priors are identical, except that tight priors were set on the periods of the signals, that is, we set a Gaussian prior as given in Table E.2 with a standard deviation 1/T_{obs} in frequency, where T_{obs} is the total observation time. The posteriors of the first and second half of the data are represented in blue and red, respectively, in Fig. D.2. The 1 sigma intervals of semi amplitude and longitude at reference time overlap in all cases, including the yearly signal.
Fig. D.2. Amplitude (K) and phase (λ_{0}) posteriors computed on the first and second half of the data for the signals. 
The coefficient of correlation with the smoothed flips sign between the first and second part of the fit, and the noise seems to exhibit slightly different time scales. In Fig. D.3, we represent the distribution of the coefficient of the smoothed and inverse timescale of the noise for the first and second halves of the data (left and right, respectively). Together with the test performed on the peak amplitudes, this suggests an evolution of the noise properties on the timescale of the data, which is potentially due to instrumental and/or stellar features.
Fig. D.3. Amplitude (K) and phase (λ_{0}) posteriors computed on the first and second half of the data for the signals. 
Appendix E: Model parameters
E.1. MCMC
To compute the uncertainties on the orbital elements, we performed a Monte Carlo Markov chain analysis. In this appendix, we define the likelihood, priors, and convergence tests used.
We assume a Gaussian likelihood. We denote the time series of radial velocities with y, the density of y knowing the parameters is of the form
Our signal model f(θ) includes Kelplerian models initialized at 2.177, 3.432, 5.198, 7.951, 12.03, 17.39, and 361 d, and a linear part, with an offset and the smoothed as defined in Appendix B.2. This one is centered and normalized by its standard deviation, so that its amplitude can be interpreted as a velocity.
The noise model includes a free white noise jitter and an exponential decay term, so that the noise model is
where η = (σ_{W}, σ_{R}, τ_{R}) are free parameters and σ_{C} is fixed to 1 m s^{−1}. In total, we have 40 parameters; θ and η have 37 and 3 components, respectively.
The prior distributions on the parameters are defined in Table E.1. For the eccentricity, we first ran the simulation with a looser prior (beta distribution with α = 1 and β = 4). Each MCMC sample was taken as an initial condition for the system. Its evolution was integrated 1 kyr in the future using the 15th order Nbody integrator IAS15 (Rein & Spiegel 2015), from the package REBOUND (Rein & Liu 2012). General relativity was included via REBOUNDx, using the model of Anderson et al. (1975). We considered the system to be unstable if any two planets have an encounter below their mutual Hill radius. From the 36 782 original samples, around 1700 passed the stability test.
Variables used for the computation of the MCMC analysis and their prior distributions.
To increase the number of effective samples and obtain reliable 3σ intervals, we reran the MCMC. We redefined the prior on eccentricity for each planet i with a beta distribution such that α_{i} = 1 and β_{i} is such that the variance of the prior is equal to the variance of the empirical distribution of the eccentricities of the points that survived the integration. We note that λ_{0} is the mean longitude at the reference epoch 57 500. For planet b, the prior T_{c} was set in accordance with TESS data.
The convergence was checked by computing the number of effective samples in each parameter chain as in Delisle et al. (2018). We find that each chain has at least 18 000 effective samples, which indicates convergence of the chain. To compute the uncertainties on values of m sin i, we took the uncertainties on the stellar mass into account and generated independent samples of Gaussian distributions with a mean and standard deviation of 1.08 and 0.1 M_{⊙}, which is in accordance with Chandler et al. (2016).
E.2. Model consistency
It might happen that the parametrized model chosen (Eq. (E.1) and (E.2)) is such that no model of this class accurately represents the data. In that case, the orbital elements are not reliable as they are computed with an incorrect model.
To check whether the model is consistent with the data, we study the residuals. Following Hara et al. (2019), we define
where y is the data, W is the inverse of the covariance matrix, f is the signal model, and are the maximum likelihood values of the parameters, as given in Table E.2. If the mode is consistent, then r_{w} should be approximately behaved as a Gaussian variable of mean 0 and variance 1.
Estimates and credible intervals of the orbital parameters for circular orbital models and a noise model including a jitter and a red noise model with exponential decay.
In Fig. E.1, we represent the normalized histogram of r_{W} (blue) and the probability density function of a normal variable (red). The histogram shows moderate asymmetry, which is inconclusive. We performed a ShapiroWilk normality test (Shapiro & Wilk 1965), yielding a pvalue of 0.1, which means that the behavior of the residuals is consistent with a normal distribution. For comparison purposes, we plotted the distribution of the non normalized residuals in Fig. E.2.
Fig. E.1. Histogram of the weighted residuals (blue) and probability density function of a normal variable (red). 
Fig. E.2. Histogram of the residuals (blue) and probability density function of a normal variable (red). 
Secondly, we searched for potential correlations in the weighted residuals. For all combinations of measurement times t_{i} > t_{j}, we represent d_{ij} := r_{W}(t_{j})−r_{W}(t_{i}) as a function of t_{j} − t_{i}. We then computed the standard deviation of the d_{ij} such that t_{j} − t_{i} is in a certain time bin. Ten such intervals are considered, with a constant length in log scale. The results are represented in Fig. E.3, where it is apparent that no significant correlations remain in the residuals. For comparison purposes, in Fig. E.4, we show the same plots for the nonweighted residuals of the model, where it appears that the dispersion of d_{ij} increase as a function of t_{j} − t_{i}.
Fig. E.3. Difference between the weighted residuals as a function of the time interval between them (blue). The red stair curves represents the standard deviation of the residuals difference in each time bin. 
Fig. E.4. Difference between the residuals (not weighted) as a function of the time interval between them (blue). The red stair curves represent the standard deviation of the residuals difference in each time bin. 
In conclusion, there is no sign of missed variance and temporal correlations in the residuals, such that the signal model of Eqs. (E.1) and (E.2) seems appropriate.
All Tables
Periods appearing in the ℓ_{1} periodogram and their false alarm probabilities with their origin, the semiamplitude of corresponding signals, and M sin i with 68.7% intervals.
Periods appearing in the ℓ_{1} periodogram of the model with the highest CV score and their false alarm probabilities.
Periods appearing in the ℓ_{1} periodogram of the model with the highest approximated evidence and their false alarm probabilities.
Variables used for the computation of the MCMC analysis and their prior distributions.
Estimates and credible intervals of the orbital parameters for circular orbital models and a noise model including a jitter and a red noise model with exponential decay.
All Figures
Fig. 1. SOPHIE radial velocity measurements of HD 158259 after outliers at BJD 2457941.5059, 2457944.4063, and 2457945.4585 have been removed. 

In the text 
Fig. 2. ℓ_{1} periodogram of the SOPHIE radial velocities of HD 158259 corrected from outliers (in blue). The periods at which the main peaks occur are represented in red. 

In the text 
Fig. 3. Radial velocity phasefolded at the periods of the signals appearing in the period analysis. From top to bottom: 2.17, 3.43, 5.19, 7.95, 12.0, and 17.4 d. 

In the text 
Fig. 4. Radial velocity phasefolded, from top to bottom, at: 360, 750, and 2000 d. 

In the text 
Fig. 5. Configuration of the system estimated from the RV at the time of conjunction estimated by TESS (BJD 58766.049072). 

In the text 
Fig. A.1. SOPHIE radial velocities and nominal 1σ error bars, the three points which were removed are shown in green. 

In the text 
Fig. A.2. measurements of HD 158259. Points in blue are conserved for the analysis; the points in red are excluded from the analysis. 

In the text 
Fig. A.3. Periodogram of the photometric data computed between on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 0.979 d (in red). The four subsequent tallest peaks are, in decreasing order, at 0.817 11.636, 1.028, and 108 d (in yellow). 

In the text 
Fig. A.4. Periodogram of the bisector span computed on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 0.985 d (in red). The four subsequent tallest peaks are, in decreasing order, at 552, 85.5, 1576, and 23.48 d (in yellow). 

In the text 
Fig. A.5. Periodogram of the computed on an equispaced frequency grid between 0 and 1.5 cycle d^{−1}. The maximum peak is attained at 2900 d (in red). Besides aliases in the one day region, there is a peak at 119 and 64 days. 

In the text 
Fig. A.6. ℓ_{1} periodogram of the APT photometric measurements. 

In the text 
Fig. A.7. ℓ_{1} periodogram of the bisector span. 

In the text 
Fig. A.8. ℓ_{1} periodogram of the . 

In the text 
Fig. B.1. ℓ_{1} periodogram corresponding to the highest approximate evidence. 

In the text 
Fig. B.2. Histogram of the values of the cross validation score. Best 20% and lowest 80% are represented in blue and orange, respectively. 

In the text 
Fig. B.3. FAPs of the peaks of the 20% best models in terms of cross validation score that have a FAP > 0.05. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 

In the text 
Fig. B.4. FAPs of the peaks of the 20% best models in terms of Laplace approximation of the evidence that have a FAP > 0.05. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 

In the text 
Fig. B.5. FAPs of the peaks of the 20% best models that have a FAP > 0.05, when adding a linear activity model fitted as a Gaussian process. The periods marked in red in Fig. 2 are represented by the blue dashed lines. 

In the text 
Fig. B.6. Raw and smoothed time series. The dark blue points represent the raw data used for the prediction. The light blue lines represent the Gaussian process prediction and its ±1σ error bars (see Eq. (B.2)). The orange points are the predicted values of the Gaussian process at the radial velocity measurement times. 

In the text 
Fig. C.1. Subsequent periodograms after a fit of a Keplerian orbit model initialized at the maximum peak of the periodogram 

In the text 
Fig. D.1. Evolution of the ℓ_{1} periodogram amplitudes of peaks at different periods. Top: planets and planet candidates. Bottom: other signals. The black dashed line indicates the transition from calibration with thoriumargon to FabryPerot. 

In the text 
Fig. D.2. Amplitude (K) and phase (λ_{0}) posteriors computed on the first and second half of the data for the signals. 

In the text 
Fig. D.3. Amplitude (K) and phase (λ_{0}) posteriors computed on the first and second half of the data for the signals. 

In the text 
Fig. E.1. Histogram of the weighted residuals (blue) and probability density function of a normal variable (red). 

In the text 
Fig. E.2. Histogram of the residuals (blue) and probability density function of a normal variable (red). 

In the text 
Fig. E.3. Difference between the weighted residuals as a function of the time interval between them (blue). The red stair curves represents the standard deviation of the residuals difference in each time bin. 

In the text 
Fig. E.4. Difference between the residuals (not weighted) as a function of the time interval between them (blue). The red stair curves represent the standard deviation of the residuals difference in each time bin. 

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.