• Top
• Abstract
• 1. Introduction
• 2. Data
• 3. Analysis of the data sets
• 4. Dynamical analysis
• 5. Conclusion
• Acknowledgments
• References
• Appendix A:
• List of tables
• List of figures

© 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, Kepler-80 (Xie 2013; Lissauer et al. 2014; Shallue & Vanderburg 2018), Kepler-223 (Borucki et al. 2011; Mills et al. 2016), and TRAPPIST-1 (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⁻¹ 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 $log{R}^{\prime }{HK}_{}$$ \log {R^\prime}_{HK} $ derived from the spectra as well as ground-based photometric data. The periodicity search is performed with a ℓ₁ 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 three-body 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 $log{R}_{hK}^{\prime }$$ \log R^\prime_{hK} $, 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⁻¹ (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.

2.2. SOPHIE radial velocities

SOPHIE is an échelle spectrograph mounted on the 193cm telescope of the Haute-Provence Observatory (Bouchy et al. 2011). Several surveys have been conducted with SOPHIE, including a moderate precision survey (3.5–7 m s⁻¹), aimed at detecting Jupiter-mass 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 M-dwarfs (e.g., Hobson et al. 2018, 2019; Díaz et al. 2019).

Since 2011, SOPHIE has been used for a survey of bright solar-type stars, with the aim of detecting Neptunes and super-Earths (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 thorium-argon lamp before barycentric Jullian date (BJD) 2458181 and a Fabry-Perot 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⁻¹. The data, which were corrected from instrumental drift and outliers (see Appendix A.1 and A.2), are shown in Fig. 1.

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 $log{R}_{hK}^{\prime }$$ \log R^\prime_{hK} $ (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 signal-to-noise 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. Ground-based photometry and ancillary indicators

If the bisector span, $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $, 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 Lomb-Scargle periodogram iteratively (Ferraz-Mello 1981; Zechmeister & Kürster 2009) as well as the ℓ₁ periodogram (Hara et al. 2017) for comparison purposes. The process is presented in detail in Appendix A.4.

The $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ periodogram presents a peak at 2900 d with a false alarm probability (FAP) of 4 × 10⁻¹². 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 ℓ₁ 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 ℓ₁ 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 cross-validation and Bayesian evidence approximations (see Appendix B). In Fig. 2, we show the ℓ₁ 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⁻¹. 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.

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⁻¹.

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 phase-folded 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.

The residuals of the nine-signals model plus white noise fit have a root mean square (RMS) of 3.1 m s⁻¹, which is higher than the nominal uncertainties of the SOPHIE data (1.2 m s⁻¹). 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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,$log{R}_{hK}^{\prime }$$ \log R^\prime_{hK} $, 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 ${95}_{-23}^{+23}°$$ {95_{-23}^{+23}}^\circ $ (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.

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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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⁻¹. We find a significantly nonzero time-scale of the noise of $5.4_{-2.6}^{+0.8}$$ 5.4_{-2.6}^{+0.8} $ 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 $7.{08}_{-0.43}^{+0.36}°$$ {7.08_{-0.43}^{+0.36}}^\circ $ (so that the transit of c cannot be detected) to $9.5_{-0.48}^{+0.57}°$$ {9.5_{-0.48}^{+0.57}}^\circ $ 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 N-body 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). 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

$$\begin{array}{c}{\psi }_{12}=3{\lambda }_2-2{\lambda }_1-{\varpi }_j,j=1,2\end{array}$$$$ \begin{aligned} \psi _{12} = 3\lambda _2 - 2\lambda _1 - \varpi _j, \;\; \; j=1,2 \end{aligned} $$(1)

would librate. Following Delisle (2017), in case of three-body resonances, the so-called Laplace angle

$$\begin{array}{c}{\varphi }_{123}=3{\lambda }_3+2{\lambda }_1-5{\lambda }_2\end{array}$$$$ \begin{aligned} \phi _{123} = 3\lambda _3 + 2\lambda _1 - 5\lambda _2 \end{aligned} $$(2)

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 time-scale 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 three-body resonances, which are seen in other resembling systems (e.g., Kepler-80, 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 $1.{09}_{-0.27}^{+0.23}$$ 1.09^{+0.23}_{-0.27} $ρ_⊕. The TESS photometry does not show other transits. There exist stable configurations of the six-planet system that are compatible with the error bars.

While many compact near-resonance 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 three-body 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 Kepler-80 (MacDonald et al. 2016). However, the latter presents three-body resonances, while HD 158259 does not. It would be interesting to investigate scenarios that can explain the differences between the two systems.

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, co-funded by CNE. G.W.H. acknowledges long-term 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. 337591-ExTrA. 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 (ANR-15-IDEX-02), through the funding of the ‘Origin of Life’ project of the Univ. Grenoble-Alpes. 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/FIS-AST/28953/2017 and PTDC/FIS-AST/32113/2017) and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalizaȩão (POCI-01-0145-FEDER-028953 and POCI-01-0145-FEDER-032113). N.A-D. acknowledges support from FONDECYT #3180063. X.Du. is grateful to the Branco Weiss Fellowship–Society in Science for continuous support.

References

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 time-series 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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 time-scale 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ and bisector measurements at the dates of the three outliers are also excluded from the analysis.

In the case of the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $, besides the removal criterion based on the MAD, the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ measured after BJD 2458000 are excluded. After this date, the values of the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ are not reliable. This is due to the change of the calibration lamp from thorium-argon to Fabry-Perot, 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $. This issue is described in greater detail in Hobson (2019), p. 92. We note that excluding some of the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $. The application of the MAD criterion on the bisector span and photometry does not lead to the removal of any point.

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 two-channel precision photometer that uses a dichroic mirror and two EMI 9124QB bi-alkali 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 $logR{\prime }_{\mathrm{HK}}^{}=-4.8$$ \log R^\prime_{\mathrm{HK}} = -4.8 $. 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).

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 ℓ₁ periodograms. The periods found are used to build candidate quasi-periodic noise models in the RV.

We first present the classical periodogram analysis. We computed the generalized Lomb-Scargle periodogram (Ferraz-Mello 1981; Zechmeister & Kürster 2009) of the ancillary indicators on a grid of frequencies spanning from 0 to 1.5 cycles day⁻¹. 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⁻¹ 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 ℓ₁ periodograms for comparison purposes.

The $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ analysis was performed after the data selection, which is presented in Appendix A.2. The periodogram (see Fig. A.5) presents a clear long-term 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⁻¹ gives signals at 3500, 119, and 32 d, with FAPs 4 × 10⁻¹², 0.18 and 0.14. The long-term signal in the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ might correspond to a magnetic cycle. In Sect. B.2, we fit a Gaussian process model to the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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⁻¹ 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 peak-to-peak amplitude is very small (a fraction of a millimagnitude) and of the same order as its uncertainty.

The ℓ₁ periodograms of the photometry, bisector span, and the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ (resp. Figs. A.6–A.8) were computed on a frequency grid from 0 to 0.95 cycles day⁻¹ to avoid the one day region. The ℓ₁ periodogram aims to express the signal as a sum of a small number of sinusoidal components. Here, the ℓ₁ periodograms present numerous peaks, which is characteristic of noisy signals, as they do not have a sparse representation in frequency. The periods appearing in ℓ₁ 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.

In conclusion, the analysis of the photometry and ancillary indicators supports the existence of a long-term magnetic cycle, appearing in the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ periodogram, with a period ≥1500 d. The period cannot be resolved due to the time-span 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

$$\begin{array}{cc}{V}_{\mathit{kl}}(\theta )& ={\delta }_{k,l}({\sigma }_k^2+{\sigma }_W^2)+{\sigma }_{\mathrm{C}}^{2}c(k,l)\\ & +{\sigma }_{\mathrm{R}}^2{\mathrm{e}}^{-\frac{|t_k-t_l|}{{\tau }_{\mathrm{R}}}}+{\sigma }_{\mathrm{QP}}^2{\mathrm{e}}^{-\frac{|t_k-t_l|}{{\tau }_{\mathrm{QP}}}}\frac{1}{2}(1+cos\left(\frac{2\pi (t_k-t_l)}{P_{\mathrm{act}}}\right))\end{array}$$$$ \begin{aligned} V_{kl}(\theta )&= \delta _{k,l} (\sigma _{k}^2 + \sigma _{W}^2) + \sigma _{\rm C}^2 c(k,l) \nonumber \\ \quad&+ \sigma _{\rm R}^2 \mathrm{e}^{-\frac{|t_k-t_l|}{\tau _{\rm R}}} + \sigma _{\rm QP}^2\mathrm{e}^{-\frac{|t_k-t_l|}{\tau _{\rm QP}}} \frac{1}{2} \left(1 + \cos \left(\frac{2\pi (t_k-t_l)}{P_\mathrm{act} } \right) \right) \end{aligned} $$(B.1)

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 quasi-periodic 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 (Foreman-Mackey 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⁻¹, τ_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 ${\sigma }_{\text{W}}^2+{\sigma }_{\text{C}}^2+{\sigma }_{\text{R}}^2+{\sigma }_{\text{QP}}^2$$ \sigma_{\rm W}^2 + \sigma_{\rm C}^2 + \sigma_{\rm R}^2 + \sigma_{\rm QP}^2 $ that is greater than the total variance of the data, and we subdivided the possible values of σ_R, σ_QP in smaller steps. The correlation time-scales 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.

For each matrix in the list of candidate noise models, the ℓ₁ 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₂₀.

Each model of CV₂₀ 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₂₀ 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₂₀.

We plotted the same Fig. B.3 and the same Table B.2 for the 20% best models in terms of evidence (E₂₀) (Fig. B.4 and Table B.3). The ℓ₁ 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₂₀, 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.

When using a LASSO-type estimator (Tibshirani 1994), such as the ℓ₁ periodogram, for a fixed dictionary (here, a fixed noise model), it is common practice to select the model with cross-validation as the solution follows the so-called 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. Long-term model

It has been found that the $log{R}_{hK}^{\prime }$$ \log R^\prime_{hK} $ has a strong long-term 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $

$$\begin{array}{c}{V}_{\mathit{kl}}(\theta )={\delta }_{k,l}{\sigma }_W^2+{\sigma }_{\mathrm{C}}^{2}c(k,l)+{\sigma }_{\mathrm{R}}^2{\mathrm{e}}^{-\frac{{(t_k-t_l)}^2}{2{\tau }_{\mathrm{R}}^2}}.\end{array}$$$$ \begin{aligned} V_{kl}(\theta ) = \delta _{k,l} \sigma _{W}^2 + \sigma _{\rm C}^2 c(k,l) + \sigma _{\rm R}^2 \mathrm{e}^{-\frac{(t_k-t_l)^2}{2\tau _{\rm R}^2}}. \end{aligned} $$(B.2)

We assume that σ_R is equal to the standard deviation of the $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ (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.

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 quasi-periodic 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⁻¹.

The planet count is initialized at 0. We compute the generalized Lomb-Scargle periodogram, and we search at which frequency its maximum is attained ω₀. We then fit a Keplerian signal initialized at ω₀, but depending on the value of ω₀, we might choose to fit ω₀ − ω_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 ω₁; we then fit two Keplerian functions initialized at ω₀ and ω₁, 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⁻¹, 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.

Overall, the periodogram results are similar to those of the ℓ₁ periodogram, the most striking difference being that in the classical approach, the 7.95 d-signal does not appear, while the 17.4 d-signal has a stronger significance. Furthermore, the 1.84/2.17 candidate does not appear. In the classical approach, the 7.95 d-signal is absorbed in the Keplerian fit initialized at 17.4 d. When performing the iterative search with sinusoids only, the 7.95 d-signal appears at the fifth iteration. The decreased significance of the 17.4 d-period in the ℓ₁ periodogram compared to the classical periodogram stems from accounting for correlated noises in the ℓ₁ periodogram. As discussed in Sect. B.3, the signals with a period greater than the time-scale 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⁻¹, 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 ℓ₁ 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 ℓ₁ 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.

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 ℓ₁-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 ℓ₁-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.

The coefficient of correlation with the smoothed $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ and inverse time-scale 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 time-scale of the data, which is potentially due to instrumental and/or stellar features.

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

$$\begin{array}{c}p(y|\theta ,\eta )=\frac{1}{{\sqrt{2\pi }}^N|V(\eta )|}{\mathrm{e}}^{\frac{1}{2}{(y-f(\theta ))}^{T}V{(\eta )}^{-1}(y-f(\theta ))}.\end{array}$$$$ \begin{aligned} p({ y}|\theta , \eta ) = \frac{1}{\sqrt{2\pi }^N | V (\eta )|} \mathrm{e} ^{\frac{1}{2} ({ y}-f(\theta )) ^T V(\eta )^{-1} ({ y}-f(\theta ))}. \end{aligned} $$(E.1)

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 $log{R}_{Hk}^{\prime }$$ \log R^\prime_{Hk} $ 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

$$\begin{array}{c}{V}_{\mathit{kl}}(\eta )={\delta }_{k,l}({\sigma }_k^2+{\sigma }_W^2)+{\sigma }_{\mathrm{C}}^{2}c(k,l)+{\sigma }_{\mathrm{R}}^2{e}^{-\frac{|t_k-t_l|}{{\tau }_{\mathrm{R}}}}\end{array}$$$$ \begin{aligned} V_{kl}(\eta ) = \delta _{k,l} (\sigma _{k}^2 + \sigma _{W}^2) + \sigma _{\rm C}^2 c(k,l) + \sigma _{\rm R}^2 {e}^{-\frac{|t_k-t_l|}{\tau _{\rm R}}} \end{aligned} $$(E.2)

where η = (σ_W, σ_R, τ_R) are free parameters and σ_C is fixed to 1 m s⁻¹. 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 N-body 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.

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 λ₀ 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

$$\begin{array}{c}{r}_w:=W{(\widehat{\eta })}^{1/2}(y-f(\widehat{\theta }))\end{array}$$$$ \begin{aligned} r_w := W(\widehat{\eta })^{1/2} ({ y} - f(\widehat{\theta })) \end{aligned} $$(E.3)

where y is the data, W is the inverse of the covariance matrix, f is the signal model, and $\widehat{\theta },\widehat{\eta }$$ \widehat{\theta},\widehat{\eta} $ 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.

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 Shapiro-Wilk normality test (Shapiro & Wilk 1965), yielding a p-value 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.

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.

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

All Figures

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.