© The Authors 2025

1 Introduction

The determination of accurate and precise masses and bulk densities of transiting small planets with radial-velocity (RV) follow-up allows us to infer their composition, under reasonable assumptions of planet interior differentiation (e.g., Zeng & Seager 2008) or mixing (e.g., Dorn et al. 2017). The ever-increasing number of well-characterized small planets has revealed a surprising variety in their densities, and hence their compositions. Generally, high bulk densities correspond to rocky planets with either terrestrial (e.g., Dressing et al. 2015) or iron-rich (e.g., Bonomo et al. 2019) interiors. Intermediate densities can be reproduced by massive icy shells (water worlds) and/or thin H/He (hydrogen/helium) atmospheres shrouding the planet rocky interiors (e.g., Mortier et al. 2018), with a strong degeneracy between the icy and gaseous mass fractions (e.g., Rogers & Seager 2010; Spiegel et al. 2014). Low densities definitively require thicker H/He atmospheres (e.g., Lissauer et al. 2011).

The fundamental goal in estimating planet compositions lies in understanding the main mechanisms that give rise to the observed compositional properties and the related well-known valley in the distribution of planetary radii around solar-type stars at R_trans ∼ 1.7–2.0 R_⊕, which likely separates prevalently rocky planets (R_p ≲ R_trans) from volatile-rich ones (R_p ≳ R_trans; Fulton et al. 2017; Van Eylen et al. 2018). Such mechanisms include planet formation at different regions of the protoplanetary disks – dry or ice-rich regions inside or beyond the water iceline, respectively (Lopez 2017; Zeng et al. 2019) – and/or at different timescales, that is, in gas-rich or gas-poor environments. Several post-formation processes that may afterwards modify the planet composition have been identified, among which (i) atmospheric photoevaporation driven by intense stellar XUV radiation, mainly in the first ∼100 Myr after formation (e.g., Owen & Wu 2013, 2017; Lopez & Fortney 2014); (ii) core-powered mass loss, that is, atmospheric escape due to the luminosity of the cooling rocky core, acting on longer timescales (∼1 Gyr; Ginzburg et al. 2018; Gupta & Schlichting 2019); and (iii) erosion due to giant impacts (e.g., Liu et al. 2015; Reinhardt et al. 2022). The observed diversity in planet compositions and the radius valley are probably caused by more than one mechanism. Nonetheless, if one of the abovementioned (formation or post-formation) mechanisms dominates over the others, it may produce a distinct slope in the position of the radius valley as a function of orbital period, which indicates a possible variation of the transition between rocky and non-rocky planets with the level of incident flux received by the planet (e.g., Lopez & Rice 2018; Martinez et al. 2019). This possible dependence of the rocky/non-rocky transition with orbital period must then be compared with the inferred compositions from bulk densities, which are expected to be rocky or volatile-rich below and above the radius valley, respectively, with the changing level of incident flux. For that purpose, more precise measurements of the bulk densities of small planets with relatively long orbital periods (hence lower insolation), P ≳ 30 d, are needed, as only a few of them have been characterized well to date.

Determining the mass and bulk density of small planets with long orbital periods can be much more difficult than for shorter-period, inner planets, for three main reasons. Firstly, the semiamplitude of the (Doppler) RV signal induced by the planet decreases as P^−1/3; secondly, wider temporal baselines of RV measurements are required to properly sample long-period signals; and thirdly, these signals may be more affected by correlated noise due to stellar magnetic activity variations and/or the presence of additional planets in multiplanet systems. A clear example of this difficulty is the planetary system K2-3, which consists of three small planets −K2-3 b, c, and d with radii R_b = 2.3, 1.8, and 1.6 R_⊕ – which orbit an early M dwarf with periods P = 10.05, 24.65, and 44.56 d, the latter in the habitable zone (Crossfield et al. 2015). Even with more than 300 HARPS-N (High Accuracy Radial velocity Planet Searcher in North hemisphere) and HARPS high-precision radial velocities, it was not possible to determine the masses/densities of the two outer planets with a precision better than 3σ (Damasso et al. 2018; see also Kosiarek et al. 2019).

Kepler-10 is a “historical” planetary system containing two transiting planets: Kepler-10 b, the first rocky planet discovered by the Kepler space mission (Batalha et al. 2011), with a radius of R_b = 1.47 ± 0.03 R_⊕ and an ultra-short orbital period (USP) of P_b = 0.837 d; and Kepler-10 c, which has both a larger radius $R_{\mathrm{b}}={2.35}_{-0.04}^{+0.09}{\mathrm{R}}_{\oplus }$ and a much longer orbital period P_c = 45.294 d (Fressin et al. 2011; Dumusque et al. 2014, hereafter D14). The planets orbit a 10.5-Gyr old GV star, which is thus expected to be magnetically quiet (average CaII activity index $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\simeq -4.97$), and is located at a distance of 186 pc.

While there has been consensus on the mass and density of the inner rocky planet Kepler-10 b, determinations of the mass and density of Kepler-10 c in the literature have been largely discrepant, illustrating that the characterization of small planets with long orbital periods may be considerably problematic, even at low levels of stellar magnetic activity. For instance, with 148 high-precision RVs collected with the HARPS-N spectrograph at the Telescopio Nazionale Galileo (Cosentino et al. 2012, 2014), D14 found a RV semiamplitude of K_c = 3.26 ± 0.36 m s⁻¹, which implies a planetary mass and density of M_c = 17.2 ± 1.9 M_⊕ and ρ_c = 7.1 ± 1.0 g cm⁻³. By using more sparse 72 RV measurements acquired with the HIRES (HIgh Resolution Echelle Spectrometer) spectrograph at the Keck telescope (Howard et al. 2010), Weiss et al. (2016, hereafter W16) found K_c = 1.09 ± 0.58 m s⁻¹, which corresponds to $M_{\text{c}}={5.7}_{-2.9}^{+3.2}{}_{\mathrm{M}}$ and ${\rho }_{\mathrm{c}}={2.4}_{-1.2}^{+1.4}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}$ (one may actually even wonder whether, with a statistical significance of 1.9σ, the signal of Kepler-10 c is truly detected in the HIRES RV data). By combining the 148 HARPS-N and the 72 HIRES RVs, W16 obtained K_c = 2.67 ± 0.34 m s⁻¹, and hence M_c = 13.98 ± 1.79 M_⊕ and ρ_c = 5.94 ± 0.75 g cm⁻³. However, this combination of HARPSN and HIRES data did not solve the issue of the evident mass discrepancy (see Fig. 5, bottom right panel in W16), but led to an average mass, which was found to be closer to the value found by D14 because the number of HARPS-N RVs was approximately twice that gathered with HIRES.

In addition to comparing the mass determinations obtained with HARPS-N and HIRES, W16 confirmed the variations in the Kepler-10 c transit times with an amplitude of ∼5 min, which were previously discovered by Kipping et al. (2015). Such transit timing variations (TTVs) indicate the presence of an additional planet in the system, Kepler-10 d, which dynamically perturbs Kepler-10 c. In an attempt to take the constraints of TTVs into account, W16 performed a preliminary combined fit of TTVs and RVs, though (i) considering for simplicity a limited number of possible orbital configurations, near 2:1 or 3:2 resonance, for Kepler-10 d, and (ii) fixing the so-called TTV “super-period” at the peak of the periodogram of the TTVs at 475.22 d (see Sect. 6 in W16). By doing so, W16 found as their best solution that Kepler-10d is in apparent 2:1 resonance with Kepler-10 c, and has a period of 101.36 d and a mass of 6.84 M_⊕.

To try and reconcile the discrepancy in the RV semiamplitude, and hence planetary mass/density of Kepler-10 c, with the HIRES and HARPS-N data, Rajpaul et al. (2017) employed a framework based on Gaussian process (GP) regression with a quasi-periodic kernel to account for possible correlated noise (e.g., Haywood et al. 2014; Grunblatt et al. 2015), by simultaneously modeling the RVs and some activity indicators, such as the $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ and bisector activity indices (Rajpaul et al. 2015). As a result, they found $K_{\mathrm{c}}={1.27}_{-0.35}^{+0.42}\mathrm{m}{\mathrm{s}}^{-1}$ and $K_{\mathrm{c}}={1.64}_{-0.34}^{+0.42}\mathrm{m}{\mathrm{s}}^{-1}$ from the HIRES and HARPS-N data, respectively, and K_c = ${1.41}_{-0.23}^{+0.25}\mathrm{m}{\mathrm{s}}^{-1}$ when combining both datasets. The latter K_c would imply $M_{\mathrm{c}}={7.37}_{-1.19}^{+1.32}{\mathrm{M}}_{\oplus }$ and ${\rho }_{\mathrm{c}}={3.14}_{-0.55}^{+0.63}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}$, which are considerably lower than those measured by both D14 and W16. From model comparison using Bayesian evidence, Rajpaul et al. (2017) also found that the three-planet model was favored over the two-planet one, and reported both the orbital period and mass of Kepler-10 d, namely P_d = 102 ± 1 d and $M_{\mathrm{d}}={5.90}_{-1.01}^{+1.70}{\mathrm{M}}_{\oplus }$ in agreement with the findings of W16. Moreover, they found a higher evidence for the GP model, and derived a period of 55.5 ± 0.8 d for the quasi-periodic variations, though it is not clear whether that period is related to the stellar rotation (Rajpaul et al. 2017).

Table 1 summarizes the RV semiamplitudes and corresponding masses of the Kepler-10 planets as retrieved in the different aforementioned works from the literature. Since it is clear from all these works that Kepler-10 is a complex planetary system, we continued to monitor it in 2017, 2019, and 2020, by practically doubling the number of HARPS-N RVs within the HARPS-N/GTO program. After extracting most of the HARPSN RVs with upgraded spectra reduction softwares (Sect. 2.2), we analyzed them with several tools (Sects. 3.3 and 3.4), also in combination with the Kepler-10 c TTVs (Sect. 3.6). In this way we improved the characterization of the Kepler-10 planetary system (Sect. 4) with regard to both its orbital architecture and the determination of precise masses, and thus possible compositions, of the Kepler-10 b and Kepler-10 c transiting planets. This in turn allowed us to put additional constraints on the formation and evolution of the Kepler-10 planetary system (Sect. 5).

2 Data

2.1 Kepler photometry

Even though the Kepler light curve was analyzed in several previous works (Batalha et al. 2011; Fogtmann-Schulz et al. 2014; D14; W16; Dai et al. 2019; Singh et al. 2022), we performed a new analysis of the Kepler data to (i) re-compute the twenty-four transit times of Kepler-10 c in short (58 s) cadence, and the first two transit epochs observed in long-cadence mode (29.4 min), which are not provided by W16; and (ii) re-derive the Kepler-10 c transit parameters.

2.2 Radial velocities

2.2.1 HARPS-N radial velocities

In total, we gathered 308 HARPS-N spectra, 55 of which before the failure of the red side of the charge-coupled device (CCD) in late September 2012 (Bonomo et al. 2014), and 253 after the replacement of the CCD. The former 55 spectra were reduced with the original HARPS-N data reduction software (DRS) version 3.7, and the latter 253 with the ESPRESSO (Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations) DRS version 2.3.5 adapted to HARPS-N (Dumusque 2021). The radial velocities were extracted from all the spectra by cross-correlating them with a G2V stellar template (e.g., Pepe et al. 2002). The new DRS pipeline allows for better long-term RV accuracy, due to a careful selection of nonsaturated thorium and argon lines that are used to calibrate the spectrograph wavelength solution. In the DRS-3.7 Kepler-10 RV data, significant variations in the flux of the original thorium-argon calibration lamp induced an artificial quadratic long-term trend of ∼5 m s⁻¹ until the replacement of the thorium-argon lamp at the beginning of June 2020, which further produced a RV offset at the level of a m s⁻¹. Those systematics disappear or are at least strongly mitigated to a level below the m s⁻¹ when using the new pipeline. In addition, the new DRS uses a novel algorithm to determine the wavelength solution that is more stable from calibration to calibration. The night-to-night RV rms offset due to wavelength solution calibration is reduced from ∼80 to ∼50 cm s⁻¹ (Dumusque 2021). However, no significant improvements from the new DRS are expected for the first 55 HARPS-N spectra taken with the old CCD, given both the relatively small amount of data and the short timespan coverage. For this reason, we used the original HARPS-N DRS (version 3.7) for those spectra.

Given the high complexity of the Kepler-10 system, as shown by the discrepancy in the Kepler-10 c mass determination from previous works, we also made use of the YARARA-v2 tool (Cretignier et al. 2021, 2022) to correct for telluric lines and subtle instrumental effects that cannot be accounted for by the new DRS pipeline, as well as for possible stellar activity variations, despite the low stellar activity level. In doing so, we had to discard 17 out of 253 HARPS-N spectra (∼7%) because they did not pass the signal-to-noise ratio (S/N) and vetting criteria of YARARA-v2.

Both the 55 DRS-3.7 and 236 YARARA-v2 HARPS-N RVs used in this work are shown in Fig. 1. They were released by Bonomo et al. (2023) along with the stellar activity indicators full width at half maximum (FWHM), contrast and bisector span of the cross-correlation function (CCF), and the CaII H & K S-index and $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$¹. The improvement achieved with YARARA-v2 can be seen in Fig. 2, which shows the increase in power of the peaks at the expected orbital periods of both Kepler-10 b and Kepler-10 c (from Kepler photometry) in the generalized Lomb-Scargle (GLS, Zechmeister & Kürster 2009) periodogram of the 236 HARPS-N RVs reduced with the new DRS+YARARA-v2 (bottom panel) compared to the new DRS (top panel). For this reason, we used the former for the 236 RVs collected from 2013 to 2021 (no remarkable improvement is instead expected for the 55 RVs gathered earlier for the same aforementioned reasons that we did not apply the new DRS for their extraction). We note that we checked for possible outliers in the HARPS-N RVs by using Chauvenet’s criterion (e.g., Bonomo et al. 2023), and found none.

2.2.2 HIRES radial velocities

The HIRES Kepler-10 spectra were recorded using different deckers that define the length of the effective slit. The first 17 HIRES observations that spanned approximately 38 days were made with the B5-decker (0.861 × 3.5 arcsec²)². The longer C2decker (0.861 × 14 arcsec²), which allows for better subtraction of night sky emission and scattered moonlight as compared to the B5-decker, was used for most of the subsequent observations.

We investigated the quality of the B5 RV set by fitting a two-Keplerian (i.e., two-planet) model to subsets of the data and analyzing the best-fit parameters. The importance of the consistency of the extracted semiamplitude in time for subsets of data is discussed in Hara et al. (2022a). For this analysis, we used PyORBIT (Malavolta 2016) by fixing the periods and transit times of planets b and c, and assuming circular orbits. The only free parameters were the semiamplitudes of the two planets and the stellar RV jitter. Since planet b has a very short orbital period, we expect to get consistent semiamplitudes for this planet across the data as long as we include a sufficient number of RVs. Fitting this model to the first 17 B5-decker observations in W16, we measured a remarkably large semiamplitude of 4.8 m s⁻¹ for planet b followed by a drop to 1.6 m s⁻¹ for the closest two-month chunk of C2-decker RVs. We also compared the B5-decker RVs with RV chunks of similar duration within the HARPS-N data with and without the YARARA-v2 correction. For this, we randomly selected two-month windows of data with at least 15 observations, to match the characteristics of the B5 RV set, and fitted the two-Keplerian model described above. The median number of included measurements was equal to 20, which is comparable to the number of B5-decker observations. The extracted semiamplitude of planet b reached around 4 m s⁻¹ for only a very few subsets of the data; none reached 4.8 m s⁻¹ as extracted for the B5-decker RVs, and we found no semiamplitude jumps as large as the one between the B5 and the C2-decker RVs. Therefore, the B5 data appear to be contaminated. To avoid the few noisy B5 measurements, which complicate a meaningful extraction of orbital information, and keep the dataset uniform, we considered only the C2 RVs.

As for the HARPS-N data, we searched for possible outliers in the HIRES C2-decker RVs with Chauvenet’s criterion, and found four measurements at the epochs 2455314.006, 2455344.978, 2456908.977, and 2456909.885 BJD_UTC, which were thus discarded.

3 Data analysis

3.1 Transit timing variations and parameters of Kepler-10 c

To compute the Kepler-10 c transit timing variations, we first fitted each transit independently by (i) using the same differential evolution Markov chain Monte Carlo (DE-MCMC; Ter Braak 2006; Eastman et al. 2013) Bayesian framework as in D14; (ii) considering intervals of the light curve centered at the predicted times of the Kepler-10 c transits from the linear ephemeris in D14, with a temporal window of twice the transit duration; (iii) adopting a two-planet transit model in case of superposition of the Kepler-10 b and Kepler-10 c transits, by fixing the parameters of the Kepler-10 b transits to the solution of D14; (iv) imposing a Gaussian prior on the stellar density, as derived from previous asteroseismic analyses of the Kepler light curve; and (v) fixing both the orbital period of planet c and the coefficients of the quadratic limb-darkening law to those found by D14 (no significant changes were noted when adopting instead Gaussian priors with the values and uncertainties given in D14). For the first two transits observed in long-cadence mode (29.4 min), we oversampled the model to 1-min samples (Kipping 2010). The Kepler-10 c transit times are given in Table 2; their variations with respect to a linear ephemeris are shown in Fig. 3, and are fully consistent with those reported by both Kipping et al. (2015) and W16.

We performed a subsequent, simultaneous DE-MCMC modeling of all the Kepler-10 c transits as in D14, after adjusting the observed-calculated (O-C) offsets from the previously computed transit times. The determined transit parameters are in excellent agreement with those from previous works, and are reported in Table 5.

3.2 Analysis of the transit timing variation signal

We analyzed the TTV signal presented in Fig. 3 to see what constraints we can obtain on the presence of non-transiting planets, and to make a comparison with the signal found in the RV analysis. A periodogram of the TTVs shown in Fig. 3 yields a peak at $P_{\text{TTV\hspace{0.17em}}}={466.2}_{-16.2}^{+17.2}$ days, which encompasses with 1σ the value given by W16. TTVs with timescales of tens to a few hundred times the orbital period of the planets are typically attributed to proximity to a mean motion resonance (MMR) with an additional planet in the system, defined by P_out/P_in = (k + q)/k, with k and q integers. Either the pair is inside the MMR, in which case the period of the TTVs scales as P(M_p/M_⋆)^−2/3 (e.g., Nesvorný & Vokrouhlický 2016), or the pair of planets is relatively close to the MMR, in which case the main TTV periodicity depends only on the orbital periods of the planets (e.g., Lithwick et al. 2012). Furthermore, short-term TTVs can be attributed to synodic “chopping” effects (e.g., Deck & Agol 2015).

First, we assumed that the perturbing planet is inside a MMR with Kepler-10 c. Given the Kepler-10 c orbital period of 45.2943 d, we estimated the mass that a perturbing planet would need to generate TTVs with a period of ∼466 days using the formulas from Nesvorný & Vokrouhlický (2016). We found that the perturbing planet should have a mass ranging between about one Saturn mass (for a 7:6 resonance, either inside or outside Kepler-10 c) and about 10 Jupiter masses (for a 2:1 resonance either inside or outside Kepler-10 c). In addition to the fact that such a massive object should have easily been found in the RV data, these configurations are also at the edge of instability (Deck et al. 2013).

We then assumed that the perturbing planet is not inside, but relatively close to, a MMR of the form P_out/P_in = (k + q)/k, with q = 1 or 2. The main TTV frequency now depends only on the orbital period of Kepler-10 c, which we fixed to 45.2943 d, and the perturbing planet, with a period called the “super period” 1/(k/P_in − (k+q)/P_out) (e.g., Lithwick et al. 2012). For each resonance considered, the perturbing planet can be either the inner or the outer planet. Then, for each of these configurations, the perturbing planet can be on either side of the MMR, since we only measure the absolute value of the super period. Supposing that P_TTV is such a super period, we can derive the orbital period of the companion, depending on the MMR it is close to. These potential companion periods are summarized in Table 3.

Interestingly, the signal at 151.06 ± 0.48 d found in the HARPS-N RVs (see Sects. 3.3, 3.4, and Table C.2) is 1σ consistent with the orbital period of ${150.50}_{-0.57}^{+0.59}\mathrm{d}$ required to generate the observed TTV period if the companion is relatively close to the 3:1 MMR with Kepler-10 c on an outer orbit. Given the relatively low uncertainty on both these periods, the probability that this is due to chance is quite low. The TTV signal therefore agrees with an orbital period of ∼151 d for Kepler10 d. In addition, the abovementioned 3:1 resonance is of second order (q = 2), and so the amplitude of the TTVs induced by the proximity to this resonance vanishes for zero eccentricities (Hadden & Lithwick 2016). This explains the relatively low TTV amplitude observed, but also implies that other TTV harmonics than the proximity to the 3:1 MMR might have comparable amplitudes.

We therefore verified if the TTV signal contains additional information on the interactions between planets c and d. Generalizing the results of Deck & Agol (2015), we know that the TTV signal expands as a function of the harmonics of the mean longitude of the perturbing planet (λ_d). We show the TTVs folded with respect to these harmonics in Fig. 4. The amplitude of each harmonic j λ_d is the sum of the contribution of the synodic (i.e., conjunction) harmonic j(λ_c − λ_d) and the effect of the MMRs of the form P_d/P_c = j/m, with m an integer. As can be seen in Fig. 4, only the third harmonic shows a significant TTV amplitude. Since it is this harmonic that contains the effect of the 3:1 MMR, we therefore conclude that the proximity to this resonance is likely dominating the TTV signal.

From this analysis, we can draw three conclusions: (i) the main periodicity observed in the TTVs of planet c can be explained by the signal at 151 d also seen in the HARPS-N RVs; (ii) since the relative proximity to this resonance generates TTVs for eccentric planets only, we expect Kepler-10 c and d to have non-zero eccentricities; and (iii) even low S/N TTVs can help confirm periodic signals found in RVs.

3.3 Periodogram analyses of radial velocities

To search for periodic signals attributable to planets other than Kepler-10 b and Kepler-10 c in the HARPS-N RVs, in view of the fact that at least one more planet is expected from the TTVs of Kepler-10 c, we performed two periodogram analyses: the classical GLS periodogram and the false inclusion probability periodogram (Hara et al. 2024).

3.3.1 Generalized Lomb-Scargle periodograms and iterative planet fitting

The peaks at the periods of Kepler-10 b and Kepler-10 c are clearly visible in the GLS periodogram of the HARPS-N RVs (Figs. A.1 and 2). To search for additional signals, we first modeled the RV signals of Kepler-10 b and Kepler-10 c as in Sect. 3.4.1, removed them from the original RV time series, and ran again the GLS periodogram on the residuals of the twoplanet model. This procedure leads to a better fit of the planet signals than the sinusoidal fit from the GLS parameters, allowing in particular a slightly eccentric Keplerian fit for Kepler-10 c. The periodogram of the residuals showed a third signal with a period of 150.8 ± 0.7 d and a theoretical false alarm probability (FAP) of 3.4 ⋅ 10⁻⁷ (see Fig. A.1). We then performed a three-planet modeling, searched for other signals with the GLS periodogram in the best-fit residuals, and found a fourth signal with a period of 83.08 ± 0.27 d and a FAP of 3.6 ⋅ 10⁻⁴ (Fig. A.1); the power at ∼26 d seen in Fig. 2 instead disappeared, being somehow related to both the signals of the three planets and the data sampling. However, the inclusion of this fourth signal was not favored by Bayesian model comparisons (Sect. 3.4), and thus it was not considered in the final modeling. None of these signals was found in the usable stellar magnetic activity indicators, such as the S-index and the bisector of the CCF.

3.3.2 False inclusion probability

Periodograms are prone to aliasing. As a consequence, a combination of planetary signals and noise can create peaks at periods that do not correspond to the planet orbital period or the stellar rotation period or one of its harmonics (Hara et al. 2017; Nava et al. 2020). To assess precisely which signals are present, it is better practice to search for several planets and an appropriate noise model simultaneously. In this section, we outline an analysis which is described in more detail in Appendix B.

We first applied the Bayesian methodology described in Hara et al. (2022b), whose end goal is to compute directly the posterior probability of the presence of a planet. By T_obs we denote the total observation timespan, and define Δ ω = 1/T_obs. We considered a grid of tightly spaced frequencies ω_k spanning from 0 to 1.5 cycles per day; for each ω_k in the grid, we computed the true inclusion probability (TIP), defined as the posterior probability that there is a planet with an orbital frequency in [ω_k − Δω, ω_k+Δ ω]. As shown in Hara et al. (2024), this detection criterion is optimal in the sense that, provided the likelihood and prior models are correct, it maximizes true detections for a given tolerance to false ones.

The computation of the TIP was done for a certain choice of priors and likelihood functions. We employed uniform priors on angles, a Beta prior on eccentricity as in Kipping (2014), and a log-uniform prior on the period. As shown in Hara et al. (2022b), the prior on the semiamplitude and noise models can lead to drastic differences in statistical significance. We used either a log-uniform prior on semiamplitude, or a mixture Gaussian model on the mass, and a white noise model. The results are shown in Fig. 5.

In all cases, we obtained very significant signals at 0.83 and 44.9 d with a false inclusion probability (FIP, equal to 1 – TIP), that is, a probability of not being present, of less than 10⁻¹². Depending on the priors, the FIP of the 150–152 d signal varied from 7 ⋅ 10⁻² to 5 ⋅ 10⁻⁶. We took a further step and computed the TIP marginalized on two noise models, white and red (see B.1 for details). To further investigate the origin of the 151 d signal and try to evaluate whether it was strictly periodic or not, we applied the methods of Hara et al. (2022a). The analysis presented in detail in Appendix B shows that the 151 d signal is compatible with a purely sinusoidal signal.

3.4 Modeling of radial velocities

After revealing planet-induced Doppler signals with periodogram analyses of the HARPS-N RV dataset of Kepler-10, we fitted them with a multi-Keplerian model by (i) varying the number of planets in the system from two to four; (ii) employing three different Bayesian modeling approaches (see Sects. 3.4.1–3.4.3) in order to evaluate the robustness of the orbital solution, given the previous discrepancies in the literature; and (iii) using a Gaussian likelihood function, which accounts for additional uncorrelated (white) noise (e.g., Gregory 2005) and/or correlated (red) noise through GP regression with quasi-periodic (QP) and/or squared exponential (SE) covariance functions (e.g., Grunblatt et al. 2015; Rajpaul et al. 2015).

The Keplerian signal of each planet was modeled with five parameters: the orbital period P; the inferior conjunction time T_c, which is equivalent to the mid-transit time for the transiting planets; $\sqrt{e}\cos (\omega )$ and $\sqrt{e}\sin (\omega )$, where e and ω are the orbital eccentricity and the argument of periastron; and the RV semiamplitude K. A zero point γ and a jitter term σ_j were also fitted for each dataset. The jitter terms were summed in quadrature to the formal RV uncertainties to take additional uncorrelated noise of unknown origin (stellar, instrumental, etc.) into account.

To check whether possible correlated noise is also present in the data and might be modeled along with the planetary signals, despite the very low stellar magnetic activity level, we further considered Gaussian Processes with the following hyperparameters: the amplitude h of the GP, the exponential decay term λ, the periodic GP term P_rot, and the inverse harmonic complexity term w. The first two hyperparameters are in common between the QP and SE covariance functions (kernels).

To be consistent among the different RV analyses, we used the same priors on the model parameters, which are summarized in Table C.1, specifically

• Gaussian priors on P and T_c of Kepler-10 b and Kepler-10 c from the modeling of the Kepler transits (D14); uniform priors on P and T_c for the signals attributable to additional non-transiting planets.

• half-Gaussian priors centered on zero with σ_e = 0.098 for the orbital eccentricities of Kepler-10 c and Kepler-10 d, following the distribution of eccentricities in multiplanet systems (Van Eylen et al. 2019). This prior avoids spurious high eccentricities, and hence dynamical instabilities due to orbit crossings; spurious eccentricities may occur when the RV semiamplitudes are comparable to the noise level as in the present case (e.g., Zakamska et al. 2011; Hara et al. 2019). For the innermost planet, Kepler-10 b, we used instead a circular orbit, because its orbit is expected to be well circularized given its extremely small (short) semimajor axis (orbital period; e.g., Matsumura et al. 2008, 2010).

• uniform priors on all the other parameters (K, γ, σ_j, and the GP hyper-parameters h, λ, P_rot and w).

We estimated the values and 1σ uncertainties of the model and derived parameters from the medians and the 15.86%–84.14% quantiles of their posterior distributions, as derived with the three different methods described below. For distributions consistent with zero, we provided only 1σ upper limits.

3.4.1 Differential evolution Markov chain Monte Carlo

The DE-MCMC method is the Markov chain Monte Carlo (MCMC) version of the differential evolution (DE) genetic algorithm (Ter Braak 2006; Eastman et al. 2013, 2019). A number of chains equal to twice the number of free parameters are run simultaneously, and perform an optimized exploration of the parameter space through the automatic choice of step scales and orientations. The sampling distribution for each step of a given chain consists in (i) randomly selecting two other chains; (ii) computing the difference of the model parameters of these two chains; (iii) determining the proposal step from (ii) and Eq. (2) in Ter Braak (2006); and (iv) accepting or rejecting the proposal step according to the Metropolis-Hastings algorithm. For the identification of the burn-in steps and the convergence and well mixing of the DE-MCMC chains, we followed the criteria proposed by Eastman et al. (2013).

We ran a DE-MCMC analysis first on the HARPS-N RVs only, and then on the combined HARPS-N and HIRES datasets. We accounted for uncorrelated (white) noise only, because the GP analysis with either QP or SE kernels did not achieve convergence of the chains, but led to very low acceptance rates and unconstrained hyperparameters, as expected from the low level of stellar magnetic activity.

The results of the DE-MCMC for the various planet models are listed in Table C. 2 for the HARPS-N RVs. We used the Bayesian information criterion (BIC; e.g., Kass & Raftery 1995; Jeffreys 1998; Mukherjee et al. 1998; Burnham & Anderson 2004; Liddle 2007) to compute the relative probabilities of the models with two, three, and four planets, and found BIC_3pl − BIC_2pl = −10.2 and BIC_4pl − BIC_3pl = +3.3. This implies that the three-planet model is favored over both the two-planet and four-planet models, being ∼165 and 5 times more likely, respectively. Moreover, an additional planet with P ∼ 83 d in the four-planet model would not be compatible with the TTVs of Kepler-10 c (see Sect. 3.6). A very important outcome is that the values of the RV semiamplitudes K of planet b and, especially, c do not vary significantly with the adopted (two- or three-planet) model, being fully consistent within 1σ (see Table C.2).

Adding the C2-decker HIRES data does not change significantly the planet parameters, as expected from their low number compared to the HARPS-N RV measurements (47 vs 291). The RV semiamplitudes agree within 1σ with the values obtained with HARPS-N only (see Table C.3), though the RV semiamplitudes of both Kepler-10 c and Kepler-10 d get slightly lower by 1σ and have slightly larger relative errors; this might be due to some HIRES instrumental effects, the investigation of which goes beyond the scope of this work. In any case, we mainly rely on the solution obtained from the modeling of the HARPS-N data, given that the correction of systematics and/or low-level activity variations with YARARA-v2 could be applied to the HARPS-N spectra only.

The best-fit with the three-planet model from the DE-MCMC analysis is shown in Fig. 6 for the HARPS-N RVs, and in Fig. C.1 for the HARPS-N and HIRES data.

3.4.2 Nested sampling through MultiNest

The Kepler-10 RVs were also modeled using the publicly available Monte Carlo nested sampler and Bayesian inference tool MultiNest v3.10 (e.g., Feroz et al. 2019), through the pyMultiNest wrapper (Buchner et al. 2014). The setup was characterized by 500 live points, a sampling efficiency of 0.3, and a Bayesian tolerance of 0.5. The Bayesian model selection was performed by comparing the values of the Bayesian evidence $\ln \mathcal{Z}$ calculated by MultiNest by assigning the same a-priori probability to each model.

We found that (i) the three-planet model is strongly favored over the two-planet model with $\ln {\mathcal{Z}}_{3pl.}-\ln {\mathcal{Z}}_{2pl.}=8.4$, and (ii) the four-planet model is weakly favored over the three-planet model, given $\ln {\mathcal{Z}}_{4pl.}-\ln {\mathcal{Z}}_{3pl.}=1.1$ and the empirical scale reported in Feroz et al. (2011); the posterior distribution of the fourth signal converged to a period of $P={82.95}_{-4.64}^{+0.20}\mathrm{d}$, encompassing the 1-yr alias of the main mode at 83 d.

3.4.3 Nested sampling through PolyChord with and without multidimensional Gaussian processes

The Kepler-10 RVs show only weak (linear) correlation with activity indicators, as expected for an old, quiet star. Specifically, the Pearson’s coefficient between RV and BIS $(\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }})$ measurements is r = 23.1%, p = 0.09 (r = −8.5%, p = 0.55) before the CCD replacement, and r =−15.0%, p = 0.021 (r = 4.0%, p = 0.54) afterwards. Nevertheless, we used the multidimensional (MD) GP framework developed by Rajpaul et al. (2015) to model RVs simultaneously with BIS and $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ observations. This GP framework assumes that any observed stellar activity signals are generated by some underlying latent function and its derivatives; this function, not observed directly, is modeled with a GP (Rasmussen & Williams 2006; Roberts et al. 2013). The framework includes components to account for convective blueshift suppression and active region evolution.

In addition to the GP component to describe correlated signals simultaneously present in RV, $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ and BIS time series, our model included possible non-interacting Keplerian terms in RVs only. Mutual orbital stability of all pairs of Keplerian orbits was checked using the criterion from Gladman (1993).

For the computation of the posterior probabilities of the model parameters, we used PolyChord, a nested sampling algorithm designed to work well even with parameter spaces with very large dimensionality (Handley et al. 2015). A short discussion illuminating its favorable properties compared with MultiNest, its direct predecessor, can be found in Hall et al. (2018); a detailed study testing PolyChord and validating its use for joint modeling of exoplanets and stellar activity in RVs is given by Ahrer et al. (2021). We called PolyChord via the pypolychord Python wrapper, leaving sampling parameters at their default values, except for the stopping (precision) criterion, which we changed from the default 10⁻³ to a (much) more stringent 10⁻⁹, to minimize the scatter in model evidence values from run to run (Ahrer et al. 2021). We ran all PolyChord sampling on a high-performance computing platform, typically using several hundred CPU cores simultaneously.

In addition to the MD GP analysis, we also used PolyChord with models that describe RVs only, without any GP component: partly as a comparison with the MD GP analysis, and partly to be able to crosscheck results obtained with other approaches (e.g., MultiNest, Sect. 3.4.2). We applied both classes of models (MD GP; no GP) to HARPS-N data only, and to HARPS-N plus HIRES data, for a total of four separate modeling ‘setups’, (i)–(iv). For each of these four setups, we explored models with between 0 and 4 Keplerian terms (each model assumed equally likely a priori), and performed three separate PolyChord runs for each model within a setup to obtain robust constraints on the evidences. The results from this analysis appear in Table C.4.

Across the four setups, planets b and c were reassuringly and resoundingly detected (Δ ln $\mathcal{Z}$ ∼ 30, compared to models without either planet). A model incorporating three Keplerians (planets b and c, plus one non-transiting planet) was strongly favored over two-planet models across all four setups, with Δ ln $\mathcal{Z}$ ∼ 6. As for the planet orbital parameters, K_b was strongly consistent while K_c was broadly consistent (within about 1σ) across all setups and applicable models. For the third Keplerian, we obtained a period of P_d ∼ 151 d (with a small secondary peak at 124.5 d) and semiamplitude K_d ∼ 1.4 m s⁻¹ across all setups and applicable models – this consistency lent additional weight to the conclusion that the third Keplerian term corresponds to a genuine planet. By contrast, there was insufficient statistical evidence for a fourth planet in any of our setups, with Δ ln $\mathcal{Z}$ ≲ 1 for the best 4-four-planet versus the best three-planet models. For a fourth Keplerian, the period posterior distributions peaked at P ∼ 83 d (corresponding semiamplitude K ∼ 0.8 m s⁻¹) with a prominent secondary peak at 114 d.

Given the strong consistency between the results we obtained with the different modeling setups, and for direct comparison with results from other methods, we present in Table C.2 the results from the MD GP analysis applied to HARPS-N data only. We note that our planet parameter posteriors were generally broader when including HIRES data, with evidence in favor of the detection of the known planets b and c weakened slightly (mirroring Rajpaul et al. 2021). We also note that the MD GP modeling with HARPS-N data only led us to a GP period of ${54}_{-9}^{+10}\mathrm{d}$, and a 1σ upper limit on the activity RV semiamplitude of <1.32 m s⁻¹. As already pointed out by Rajpaul et al. (2017), it is nevertheless unclear whether the GP period of ∼54 d corresponds to the stellar rotation period, given that neither the Kepler light curve nor the time series of the CCF and CaII activity indicators show clear signals attributable to the stellar rotation.

3.5 Search for Kepler-10 d transits

Given the presence of a significant signal at ∼151 d, we analyzed the Kepler light curve to check whether any transit could have been missed in previous studies of this system. In particular, we first modeled the light curves using the python package wotan (Hippke et al. 2019) and then performed a transit least squares (TLS) periodogram (Hippke & Heller 2019) on the residuals of the light curve, after modeling the transits of both Kepler-10 b and Kepler-10 c. We found no significant periodicity at ∼151 d or elsewhere, despite the Kepler curve being about ∼1480 d long (minus a few gaps in between). We also visually inspected the light curve around the eight predicted transits of Kepler-10 d (from the ephemeris of Table C.2) and found no significant dip. We point out that, based on the estimated RV semiamplitude of the outer planet, we would expect transits of a comparable depth to planet c (∼0.5 mmag), and thus easily detectable (see Fig. 7). Planet d would transit if it had the same inclination as planet c, since b_d ≤ 1: $$b_{\mathrm{d}}=\frac{a_{\mathrm{d}}}{R_{\star }}\cos i_{\mathrm{d}}\left(\frac{1-e_{\mathrm{d}}^2}{1+e_{\mathrm{d}}\sin {\omega }_{\mathrm{d}}}\right)={0.64}_{-0.08}^{+0.07}.$$(1)

However, the transiting condition breaks with a small mutual inclination of only 0.2 deg between the two planets.

3.6 Simultaneous modeling of radial velocities and transit timing variations

The TTVs of planet c have an amplitude, A_TTV³, of about 5.6 minutes and show a quasi-sinusoidal pattern (see Fig. 3), which suggests an interaction with an additional non-transiting planet on an external orbit. For this reason, we decided to run a dynamical analysis with the latest version of TRADES⁴ (Borsato et al. 2014, 2019; Nascimbeni et al. 2023) and simultaneously fit the transit times (T_c s) and the RVs during the numerical integration of the planetary orbits. The computational time required by the dynamical analysis scales with the number of orbits of the inner planet with respect to the integration time (about 4222 days). This implies that the inclusion of planet b (P_b = 0.837 d) causes the computational time to increase disproportionately. Therefore, we had to remove its signal from the RV dataset, according to the solution from Sect. 3.4.1. In addition to the considerable reduction in computation time, this choice is justified by the fact that the USP planet b cannot induce a TTV on planet c, so we assumed a three-body system with the star, planet c, and the hypothetical planet d.

The TRADES code cannot manage complex stellar activity, it can only fit an RV offset (γ) and a jitter term (log₂ σ_j) for each RV dataset; this is not an issue in our case, given the low level of stellar activity and the fact that the planet parameters do not critically depend on the noise model (Sect. 3.4.3). We fixed the inclination of planet c (i_c) to 89.6 deg (Table 5), and assumed the external planet d to have an inclination of 90 deg. We also fixed the longitude of nodes (Ω) of both planets to 0 deg. We then varied the ratio of the planetary to stellar mass (M_p/M_⋆), the period (P), the eccentricity (e) and the argument of periastron (ω) as $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$, and the mean longitude (λ⁵), of both planets c and d. All the orbital parameters are osculating astrocentric parameters at the reference time BJD_TDB−2450000.0 = 7081.0 and have uniform priors with physical boundaries, which are: M_{c, d} = $\mathcal{U}(0.1,20)$ M_⊕, P_c = $\mathcal{U}(42,48)$ d, P_d = $\mathcal{U}(50,250)$ d, e_{c, d} = $\mathcal{U}(0.0,0.5)$, ω_{c, d} = $\mathcal{U}(0.0,360)$ deg, λ_{c, d} = $\mathcal{U}(0.0,360)$ deg. Firstly, we searched for an optimal configuration with a differential evolution algorithm (DE, Storn & Price 1997) implemented as PyDE within pytransit⁶ (Parviainen 2015) with a population of 76 configurations that evolved for 76000 generations. Then, we passed the optimal configuration to the python code emcee (Foreman-Mackey et al. 2013, 2019). We used as walkers the same number of the DE population (76), and we mixed the sampler algorithm as in Nascimbeni et al. (2023). We ran emcee for 2.5 million steps, removed as burn-in the first two million, and used a thinning factor of 100. We computed as parameter uncertainties the high density interval (HDI) at 68.27%⁷ of the posterior distribution, and our best-fit solution as the maximum likelihood estimator (MLE), that is, the configuration of parameters that maximizes the log-likelihood (log $\mathcal{L}$) of the posterior distribution within the HDI. Table 4 gives the results of the best-fit parameters, and Figs. 8 and 9 show the best-fit models of TTVs and RVs, respectively. We checked the best-fit model assuming that planet d does not transit, and found an average difference on the transit times of less than one second. This is well below our sensitivity threshold given by the mean transit uncertainty of two minutes. Therefore, the effect of using an orbital inclination slightly different from 90 deg for planet d is negligible in our analysis.

The parameters are in agreement with the different solutions of Sect. 3. In particular, the eccentricities of both planets c and d are more precise than those retrieved in the DE-MCMC and MultiNest analyses (Sects. 3.4.1 and 3.4.2). We decided to test a model with the same number of bodies, but fixing e_d = 0, and we ran PyDE and emcee with the same number of walkers (76), but with a lower number of steps (76 000 PyDE generations, one million emcee steps with 400 000 as burn-in, and a further thinning factor of 100). We computed the BIC for the eccentric (ecc) and for the circular (circ) best-fit model, and we found that the eccentric model is strongly favored: Δ BIC = BIC(ecc) − BIC(circ) = − 13. (e.g., Kass & Raftery 1995).

4 Planet parameters

We updated the orbital and physical parameters of the Kepler-10 planets in Table 5, based on (i) the stellar parameters in D14; (ii) the results of the DE-MCMC transit and RV modeling (Sect. 3), also for consistency with both D14 and Bonomo et al. (2023) (we recall that fully compatible results were obtained with the other two techniques employed; see Table C.2); and (iii) the combined analysis of TTVs and RVs with TRADES for planets c and d.

Thanks to the collection, refined extraction (Sect. 2.2), and analysis (Sect. 3) of nearly 300 HARPS-N RVs, we determined the masses of Kepler-10 b and Kepler-10 c to be M_b = 3.24 ± 0.32 M_⊕(10 σ precision) and M_c = 11.29 ± 1.24 M_⊕ (9 σ precision), respectively (see Table 5, DE-MCMC solution), the latter being approximately in the middle between the previously determined values of M_c ∼ 6–7 M_⊕ (W16, Rajpaul et al. 2017) and M_c ∼ 17 M_⊕ (D14). The refined corresponding densities ρ_b = 5.54 ± 0.64 g cm⁻³ and ρ_c = 4.75 ± 0.53 g cm⁻³ allow us to better infer the planet compositions (Sect. 5), and place both Kepler-10 b and Kepler-10 c among the best characterized small planets.

For the outer non-transiting planet Kepler-10 d, we found a more accurate orbital period of P_d = 151.06 ± 0.48 d from the HARPS-N RVs (Table 5), which is consistent with both the analysis of TTVs (Sect. 3.2) and the combined analysis of TTVs and RVs (cf. Sect. 3.6 and Table 5), and derived a minimum mass of M_d sin i = 12.00 ± 2.15 M_⊕.

5 Discussion and conclusions

Figure 10 shows the position of both Kepler-10 b and Kepler-10 c in the radius-mass diagram of small exoplanets with mass and radius determinations better than 4σ and 10σ, respectively.

The innermost planet Kepler-10 b is a rocky super-Earth that likely formed inside the water snowline at ∼1−3 AU, that is, in a possibly dry environment, and then migrated inward to the current semimajor axis where it is observed now. The planet may possibly have lost a primordial atmospheric H/He envelope because of the strong XUV stellar irradiation received in the first ∼100 Myr (e.g., Owen & Wu 2013), but may have retained a thin, recently formed collisional secondary atmosphere (Singh et al. 2022). Since Kepler-10 b is located on the pure-silicate isocomposition curve (Fig. 10, left panel), its core mass fraction is broadly consistent with zero, though a small iron core may be present.

On the contrary, Kepler-10 c, by its position in Fig. 10 and its low equilibrium temperature (T_eq ≃ 580 K), may be a water world, with a significant (∼40–70%) mass fraction of H₂O in a differentiated or miscible interior (e.g., Luo et al. 2024). Unlike many of the exoplanets in the figure, which are above the pure H₂O mass-radius curve due to the presence of some atmospheric envelope, Kepler-10 c lies below it, despite being on the right-hand side of the Cosmic Hydrogen and Ice Loss Line, where atmospheric loss due to stellar radiation is expected to be negligible (cf. Zeng & Jacobsen 2024). The existence of water worlds is predicted by models of planet formation and migration (e.g., Venturini et al. 2024; Burn et al. 2024). Kepler-10 c may thus join the fifteen or so possible water worlds that orbit relatively far from their FGK dwarf hosts, and hence have rather low equilibrium temperatures T_eq ≤ 600 K, located in between the 100% silicate and 100% water isocomposition curves (Fig. 10, right panel; see also Luque & Pallé 2022 for possible water worlds around M dwarfs).

Nevertheless, due to the well-known degeneracy in planet compositions from the measurement of the bulk density, a very thin H/He atmospheric envelope of less than 1% of its total mass on top of a rocky interior (instead of a water-world composition) cannot be excluded. In this case, Kepler-10 c might have accreted a more massive H/He envelope and lost a large fraction of it through “boil-off” (Owen & Wu 2016) and/or, more likely, core-powered mass loss (Gupta & Schlichting 2019).

The rather large density of Kepler-10 c is in agreement with the trend that sub-Neptunes that are relatively far from MMR⁸ are denser (Leleu et al. 2024). This may provide an additional potential scenario to explain its density: the ejection of all or most of its atmosphere in post-disk dynamical instabilities (e.g., Izidoro et al. 2017).

Kepler-10 d might be similar to Kepler-10 c in terms of its composition, its minimum mass being consistent with the mass of Kepler-10 c within 1σ. Kepler-10 d might thus potentially be a water world or have a composition more similar to Uranus and Neptune, that is, with a larger H/He envelope, given that it is colder and likely formed farther away. In any case, the lack of the measurement of both its radius and true mass makes any discussion quite speculative.

Given the long RV time series of Kepler-10 spanning ∼11 years, we estimated the sensitivity of our data (completeness) to the presence of possible outer cold Jupiters as in Bonomo et al. (2023, see their Sect. 3.3). The completeness is shown in Fig. 11, where yellow cells indicate 100% detectability of the cold Jupiter, lighter to green cells show decreasing detection probability, and the blue ones correspond to 0% detectability. It is clear that there are no giant planets orbiting the Kepler-10 star within 10 AU, as is the case for the great majority of Kepler and K2 small planet systems followed by HARPS-N (Bonomo et al. 2023); otherwise, we would have easily detected it. This is consistent with the scenario that no giant planet could act as a dynamical barrier to the inward migration of the sub-Neptunes Kepler-10 c and Kepler-10 d (e.g., Izidoro et al. 2015).

This work highlights the crucial importance of intensive long-term RV monitoring of planetary systems with relatively long-period (P ≳ 30 d) transiting small planets, also in view of the forthcoming PLAnetary Transits and Oscillations of stars (PLATO) space mission (Rauer et al. 2014, 2024). As in the case of Kepler-10, even several hundreds of RV measurements may be needed to (i) precisely determine the masses and hence bulk densities of the transiting planets, and (ii) properly characterize the system architecture by searching for additional non-transiting small and giant planets in outer regions. The increasing number of planetary systems better characterized in this way will lead to a more comprehensive understanding of the formation and evolution of sub-Neptunes and super-Earths.

Acknowledgements

The authors would like to thank the anonymous referee for her/his comments, which allowed them to improve the manuscript. They also thank Lauren Weiss for helpful discussions about the HIRES RVs and the deckers used. This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias (GTO programme). The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard-University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), the University of St. Andrews, Queen’s University Belfast and the University of Edinburgh. This paper is based on data collected by the Kepler mission. Funding for the Kepler mission was provided by the NASA Science Mission directorate. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. The authors acknowledge financial contribution from the INAF Large Grant 2023 “EXODEMO” and from the SNSF. L.Bo. acknowledges support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0. L.Z. acknowledges the support from the DOE-NNSA grant DE-NA0004084 awarded to Harvard University: Z Fundamental Science Program (PI: Stein B. Jacobsen). M.C. acknowledges the SNSF support under grant P500PT_211024. A.L. acknowledges support of the Swiss National Science Foundation under grant number TMSGI2_211697. X.D. acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No. 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No. 200021_215200). A.M. acknowledges funding from a UKRI Future Leader Fellowship, grant number MR/X033244/1 and a UK Science and Technology Facilities Council (STFC) small grant ST/Y002334/1. R.D.H. is funded by the UK Science and Technology Facilities Council (STFC)’s Ernest Rutherford Fellowship (grant number ST/V004735/1). This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch.

Appendix A Generalized Lomb-Scargle periodograms of the HARPS-N radial velocities

Appendix B Exploring the parameter space

B.1 False inclusion probability

This appendix gives the details of the analysis of the HARPSN data on Kepler-10 processed with YARARA-v2 (Cretignier et al. 2021, 2022). We define the true inclusion probability (TIP) as the posterior probability of the following event: for a given range of periods I, there is at least one planet with period P ∈ I. By analogy with the FAP, we also define the false inclusion probability (FIP), which is the probability that there is no planet with period P in interval I. Formally, denoting the data by y, the TIP and FIP are $${\mathrm{T}\mathrm{I}\mathrm{P}}_I=\mathrm{Pr}\{\mathrm{\exists }P,P\in I\mid y\}.$$(B.1) $${\mathrm{F}\mathrm{I}\mathrm{P}}_I=1-{\mathrm{T}\mathrm{I}\mathrm{P}}_I.$$(B.2)

In practice, we can evaluate Eq. (B.2) as follows. We suppose that there is a maximum number of Keplerian signals in the data k_max. We denote by p(θ∣k) the prior probability of model parameters knowing there are k planets, and p(y∣θ, k) the likelihood of the data y knowing the number of planets k and the orbital parameters. We suppose a prior probability for the model with k planets Pr{k}. Then, as defined in Eq. (B.1), the TIP is $${\mathrm{TIP}}_I={\displaystyle \sum _{k=0}^{k_{\text{max\hspace{0.17em}}}}}\mathrm{Pr}\{\mathrm{\exists }i\in [1..k],P_i\in I\mid y,k\}\mathrm{Pr}\{k\mid y\},$$(B.3)

where the Pr{k∣y} term is computed from the Bayesian evidences of k planet models, and is called PNP for posterior number of planets, as defined in Brewer (2014). It is the probability to have k planets knowing the data. Denoting by ${\mathcal{Z}}_k$ the Bayesian evidence of the model with k planets, assuming there are at most k_max planets, and denoting by Pr{i} the prior probability to have i planets, $$\mathrm{Pr}\{k\mid y\}=\frac{{\mathcal{Z}}_k\mathrm{Pr}\{k\}}{{\displaystyle \sum _{i=0}^{k_{max}}}{\mathcal{Z}}_i\mathrm{Pr}\{i\}}.$$(B.4)

The main advantage of the TIP/FIP is that it is very easy to interpret. If the model is correct, among detections made with TIP = α = 1 – FIP, on average a fraction α are correct detections and a fraction 1−α are spurious ones.

We perform the FIP calculation up to four planets with two different types of priors for the semiamplitude. The rationale is that, as we have shown in Hara et al. (2022b), the significance of low-amplitude signals highly depends on the chosen prior. In the two sets of priors, we use a log-uniform prior on period, we constrain the periods and phase of two of the planets according to the constraints from the transit. The eccentricity has a Beta prior as in Kipping (2014), the prior on the argument of periastron is uniform.

We then compute the FIP/TIP as the probability to have at least one planet with frequency in the interval [f−1/T_obs, f + 1/T_obs] where T_obs is the total time of observation, for an array of equispaced frequencies. In the first set of priors, we use a loguniform prior from 0.1 to 20 m s⁻¹ on the RV semiamplitudes, and obtain the FIP periodogram shown in Fig. 5(a).

In the second set of priors, the prior is not defined on semiamplitude but on the mass. This means that, for each set of period, eccentricity and time of passage at periastron (or equivalently, mean anomaly at a reference time), the prior on semiamplitude is scaled to correspond to the same prior on mass. This choice implicitly supposes that the probability of occurrence of a given type of planet (super-Earth, Neptune) does not depend on semimajor axis. Since, for a given mass, the semiamplitude decreases with separation, this will also have the effect to boost the significance of low amplitude signals. We use a Gaussian mixture model for the prior on mass, which follows a Rayleigh prior with σ = M_⊕ or σ = 15 M_⊕ with a probability 1/2. With this choice, we can marginalize the likelihood analytically over the linear parameters as shown in Appendix C of Hara et al. (2022b). In that case, a 68 d signal reaches a FIP of 4% (see Fig. 5, b).

We did the same analyses as before but including a red noise term with an exponential kernel (free jitter and time-scale). The corresponding FIP periodograms are shown in Fig. B.1, where the 150.7 d signal is much less significant and completely disappears, respectively. However, we caution that the red noise model is expected to damp the amplitude of long-period signals, as shown in Delisle et al. (2020). It is common for this type of noise model to suppress viable planets. If what happens at low frequency was truly red noise, then we would expect the white noise model to find several unclear planet candidates at long periods, due to a diffuse power at low frequencies. However, in Kepler-10 the white noise model shows a clear candidate with a well defined period. As a consequence, the interpretation of the FIP periodogram with a red noise model with exponential kernel is subject to caution.

In Table B.1 and Table B.2 we show the evidences $\mathcal{Z}$ and uncertainty on them (standard deviation of three runs) for the log-uniform prior on K, white and red noise assumptions, respectively. We can now compute the TIP of the 150.7 d planet averaged over the noise model. Denoting by I the frequency interval $[\frac{2\pi }{150.7}-\mathrm{\Delta }\omega ,\frac{2\pi }{150.7}+\mathrm{\Delta }\omega ]$, we have $$p(\omega \in I\mid y)=p(\omega \in I\mid W,y)p(W\mid y)+p(\omega \in I\mid R,y)p(R\mid y),$$(B.5)

where W and R represent the white and red noise signal. The two terms p(ω ∈ I ∣ W, y) and p(ω ∈ I ∣ R, y) are equal to 1-FIP of the 150.7 d signal for the white and red noise models, respectively. We have p(ω ∈ I ∣ R, y) = 0. Assuming white and red noise have both a prior probability of 1/2, we have $$p(W\mid y)=\frac{p(y\mid W)}{p(y\mid W)+p(y\mid R)}$$(B.6)

and $$p(y\mid W)={\displaystyle \sum _{k=0}^3}p(y\mid W,k)p(k)$$(B.7)

where k is the number of planets. With the values of Table B.1 and Table B.2, assuming p(k) = 1/4, we have $$p(y\mid W)={\mathrm{e}}^{-732.70}/4$$(B.8) $$p(y\mid R)=({\mathrm{e}}^{-734.77}+{\mathrm{e}}^{-735.53})/4$$(B.9)

which means that p(W ∣ y) = 0.84. As a result, the probability to have a signal at 150.7 d is (1–10^−2.16) × 0.84 ∼ 0.84. The probability that there is a signal in I marginalized on the noise models considered here is 84%. Conversely, we have a FIP for the 150.7 d signal of 1–0.84 = 0.16.

B.2 Apodized sine periodograms, white noise

The result from the two previous analyses show that there is a significant signal at 150.7 d, and hints of signals at 24 and 68 d, the latter being the most promising, but still not significant enough to be a clear detection. The question is then mainly to determine if the 150.7 d signal, which is clearly significant, can be attributed to a planet.

To address this question, we analyze the data in a different framework: the apodized sine periodogram (Hara et al. 2022a). This framework, proposed by Gregory (2016), is designed to determine if signals are present throughout the whole dataset, and have a constant phase, amplitude and frequency. This approach is particularly useful when there is no clear model for the data. In practice, we fit the apodized sine function $${\mu }_K(t,\omega ,\tau ,t_0,A,B)=w(\tau ,t_0)(A\cos \omega t+B\sin \omega t),$$(B.10)

where w(τ, t₀) is the apodization function, t₀ is its center and τ its temporal width. The resulting signal μ is then wavelet-like, quasi periodic from ∼t₀−τ to ∼t₀+τ and 0 elsewhere. We here use Gaussian and box-shaped functions, that is $$w_G(\tau ,t_0):={\mathrm{e}}^{-\frac{{(t-{\tau }_0)}^2}{2{\tau }^2}}$$(B.11) $$w_B(\tau ,t_0):={\mathbf{1}}_{[t_0-\frac{\tau }{2},t_0+\frac{\tau }{2}]}(\tau ,t_0)$$(B.12)

We then compute a periodogram, using the same defintion as Baluev (2008), $$z(\omega ,t_0,\tau )={\chi }_H^2-{\chi }_K^2(\omega ,t_0,\tau ),$$(B.13)

where ${\chi }_H^2$ is the χ² of the residuals after fitting the model H (the base model, without sinusoidal function) and the model K(t, ω, τ, t₀), that is, linearly fitting A and B for model Eq. B. 10 simultaneously with the base model. Then, for a given value of τ, we represent z as a function of ω maximized on t₀, that is, for a given frequency and timescale, the best fitting model when the center of the window varies.

The transiting planets are added to the base model (they are assumed to be there in the null hypothesis). We take four values of τ in the grid, such that τ/T_obs = 10, 1/3, 1/9, 1/27 where T_obs is the total observation time-span. The apodized sine periodogram for Kepler-10 (with the transiting planets in the null hypothesis model) is shown in Fig. B.2. The best fit occurs at τ = T_obs/3, and the corresponding model is shown in Fig. B.3(b). Given the gap in the data, the best fitting timescale might be spuriously shorter than the true ones due to random fluctuations. The right panel of Fig. B.2 shows a statistical test to see if some values of τ can be excluded (see Hara et al. 2022a). It shows that the highest value of τ is still 2σ compatible with the data.

In Fig. B.3(a) we show the evolution of the phase and amplitude of that signal over the observations. While the phase is constant, the amplitude tends to vanish towards the end of the observations. We then multiplied the sinusoidal model of the signal by a so called apodization factor (Gregory 2016; Hara et al. 2022a), e^{(t−t₀)2/(2τ2)} where t is the time, and t₀ and τ are free parameters. The best fit is represented in Fig. B.3(b), where it seems that the amplitude of the signal vanishes at the end of the dataset. As said above, this is simply the best fitting timescale, but a purely sinusoidal signal (τ = 10 T_obs) is still compatible with the data. We performed the same analysis on the ancillary indicators and did not find trace of signals at ∼151 d.

Appendix C Radial-velocity analyses

References

---

All Tables

All Figures