# Architecture of planetary systems with and without outer giant planets

J.-B. Delisle, J P Faria, D. Ségransan, E. Fontanet, W. Ceva, D. Barbato, S G Sousa, N. Unger, A. Leleu, F. Bouchy, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

J.-B. Delisle, J P Faria, D. Ségransan, E. Fontanet, W. Ceva, et al. Architecture of planetary systems with and without outer giant planets. Astronomy & Astrophysics - A&A, 2025, 703, pp.A78. 10.1051/0004-6361/202555445. [hal-05470458](https://hal.science/hal-05470458v1)

HAL Id: hal-05470458

https://hal.science/hal-05470458v1

Submitted on 21 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Architecture of planetary systems with and without outer giant planets

## I. Inner planet detections around HD 23079, HD 196067, and HD 86226

J.-B. Delisle<fnref n="1" />,* D. J. P. Faria<fnref n="1" />, D. Ségransan<fnref n="1" />, E. Fontanet<fnref n="1" />, W. Ceva<fnref n="1" />, D. Barbato<fnref n="2" />, S. G. Sousa<fnref n="11" />,<fnref n="12" />, N. Unger<fnref n="1" />, A. Leleu<fnref n="1" />, F. Bouchy<fnref n="1" />, M. Cretignier<fnref n="3" />, R. F. Díaz<fnref n="4" />,<fnref n="5" />, X. Dumusque<fnref n="1" />, Y. G. C. Frensch<fnref n="1" />, N. C. Hara<fnref n="6" />, G. Laughlin<fnref n="7" />, G. Lo Curto<fnref n="8" />, C. Lovis<fnref n="1" />, M. Marmier<fnref n="1" />, M. Mayor<fnref n="1" />, L. Mignon<fnref n="1" />,<fnref n="9" />, C. Mordasini<fnref n="10" />, F. Pepe<fnref n="1" />, N. C. Santos<fnref n="11" />,<fnref n="12" />, and S. Udry<fnref n="1" />

Received 8 May 2025 / Accepted 5 September 2025

## ABSTRACT

Understanding the link between outer giant planets (OGPs) and inner light planets (ILPs) is key to understanding planetary system formation and architecture. The correlation between these two populations of planets is debated both theoretically – different formation models predict either a correlation or an anticorrelation – and observationally. Several recent attempts to constrain this correlation have yielded contradictory results, due to small-number statistics and heterogeneous samples. We present an ongoing long-term observational effort with CORALIE, HARPS, and ESPRESSO to probe the ILP occurrence in systems with and without OGP. In this first article of a series, we discuss how, from the design to the observations, we ensured the homogeneity of the samples, both in terms of stellar properties and observing strategy. We also present the first three detections of ILPs in our OGP host sample. We find a 8.3 $m_⊕$ planet at 5.75 d around HD 23079, a 10.4 $m_⊕$ planet at 4.6 d around HD 196067, and we confirm the 7.5 $m_⊕$ planet at 3.98 d around HD 86226. While a rigorous statistical analysis of our samples will be performed in subsequent studies, the relatively low number of detections in our sample seems to contradict previous studies that found a strong OGP-ILP correlation.

**Key words.** techniques: radial velocities – planets and satellites: general

## 1. Introduction

Characterizing the global architecture of planetary systems is crucial for understanding their formation history and dynamical evolution. Although more than 5000 exoplanets have already been detected<fnref n="1" />, it remains extremely challenging to fully characterize both the inner and outer regions of planetary systems due to the disparate planetary orbital periods, ranging from hours to decades, and the detection limits and biases of the different observing methods (transit, radial velocities, imaging, etc.). Thus, the coupling between the inner and outer parts of planetary systems during planet formation and subsequent evolution remains poorly constrained by observations. In particular, the

question of whether the presence of outer giant planets (hereafter OGPs) enhances or inhibits the formation of inner light planets (hereafter ILPs) remains open, both theoretically and observationally.

The main planet formation models have different theoretical predictions. In the inward migration theory, in which ILPs are initially formed in the outer part of the protoplanetary disk and migrate to the inner part due to disk-planet interactions, the presence of an OGP is proposed to hinder inward migration (e.g., Izidoro et al. 2015). Similarly, in the pebble accretion theory, ILPs form in the inner part of the disk owing to a flux of pebbles drifting inward from the outer region. The presence of an OGP would, in that case, inhibit the pebble flux and thus the formation of ILPs (e.g., Rice et al. 2006; Ormel et al. 2017; Lambrechts et al. 2019). Thus, these formation models predict an

<footnote n="1">¹ Département d'astronomie, Université de Genève, chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="2">² INAF - Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="3">³ Department of Physics, University of Oxford, Oxford OX13RH, UK</footnote>

<footnote n="4">⁴ Instituto Tecnológico de Buenos Aires (ITBA), Iguazú 341, Buenos Aires, CABA C1437, Argentina</footnote>

<footnote n="5">⁵ Instituto de Ciencias Físicas (CONICET / ECyT-UNSAM), Campus Miguelete, 25 de Mayo y Francia, (1650) Buenos Aires, Argentina</footnote>

<footnote n="6">⁶ Aix Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="7">⁷ Department of Astronomy, Yale University, 52 Hillhouse Avenue, New Haven, CT 06511, USA</footnote>

<footnote n="8">⁸ European Southern Observatory, Av. Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile</footnote>

<footnote n="9">⁹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="10">¹⁰ Division of Space Research and Planetary Sciences, Physics Institute, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="11">¹¹ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="12">¹² Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote>* Corresponding author: jean-baptiste.delisle@unige.ch
¹ (e.g., https://exoplanet.eu).</footnote>
anticorrelation between ILP and OGP populations. However, in the in situ formation model, orbital migration is assumed to play a limited role, and the protoplanetary disk is initially assumed to be sufficiently massive for ILPs to form in the inner region by accretion of solids and gas in their vicinity. Such a massive disk also tends to form gas giant planets in its outer region (e.g., Chiang & Laughlin 2013; Schlecker et al. 2021). Thus, a correlation between ILPs and OGPs would be expected in the in situ formation model.

Several recent studies have attempted to constrain this correlation from observations, leading to contradictory results. Most of these studies rely on comparing the occurrence rates of OGPs in blind radial velocity (RV) surveys and in RV follow-up of transiting ILP hosts. The first attempts (Zhu & Wu 2018; Bryan et al. 2019) find an enhanced probability (by a factor of three, albeit with large uncertainty) of OGP around transiting ILP hosts compared to the general population (OGPs in blind RV surveys). From this, they conclude that stars hosting an OGP have about a 90% chance of also hosting ILPs.

At the same time, a HARPS (High Accuracy Radial velocity Planet Searcher) search for inner planets in 20 OGP hosts (Barbato et al. 2018) yielded no ILP detection. However, the relatively small sample of stars and the low number of HARPS measurements per star (27 on average) did not allow the authors to draw firm conclusions. A statistical analysis of RV archival data from the California Legacy Survey (CLS, Rosenthal et al. 2021), shows a nonsignificant (within $1-\sigma$) correlation between OGPs and ILPs (Rosenthal et al. 2022). These approaches (Barbato et al. 2018; Rosenthal et al. 2021), as well as the one proposed here, are direct measurements of $P$(ILP|OGP), in contrast to the indirect method proposed by Zhu & Wu (2018); Bryan et al. (2019) which relies on Bayes' theorem to derive $P$(ILP|OGP) from $P$(OGP|ILP).

More recently, Bonomo et al. (2023) present a dedicated HARPS-N follow-up of 38 transiting ILP hosts to search for OGPs. They report a lower but compatible rate of OGP than in the general population, in contradiction with Zhu & Wu (2018) and Bryan et al. (2019). They attribute this discrepancy to the possible confusion between OGPs and magnetic cycles or binary stars. Another possible reason for this discrepancy, proposed by Zhu (2024), is that the stars observed by Bonomo et al. (2023) are metal-poor compared to the general population. However, the impact of stellar metallicity on the (anti-)correlation between the two populations is still debated (Bryan & Lee 2024; Van Zandt & Petigura 2024; Bryan & Lee 2025; Van Zandt et al. 2025; Bonomo et al. 2025).

More generally, the question of whether these populations correlate or anticorrelate is difficult to settle using archival data. Indeed, all of these studies suffer from small-number statistics (small stellar sample sizes) and/or nonhomogeneous samples, especially when comparing results from blind RV surveys and transit follow-ups. As shown by Van Zandt & Petigura (2024), some results strongly depend on the exact definitions of ILP and OGP (in terms of mass, radius, periods, and semi-major axis). This may partly reflect true physical mechanisms, but could also arise from small-number statistics and biased samples.

Additionally, a subtler issue, which we term “observer-excitement” bias, can affect some of these results. In many blind RV surveys, targets showing promising signals tend to be monitored with higher cadence, while others receive fewer observations. Consequently, blind surveys typically have more RV measurements (sometimes orders of magnitude more) for targets with planets than for those without. This bias can have two relevant effects on OGP-ILP correlations studies:

- Low-amplitude planets are more likely to be detected in multi-planetary systems where an easier-to-detect planet has already been found, increasing the likelihood of intensive monitoring;

- When correcting for completeness, systems with planets typically have more points, and thus lower detection limits, giving them greater weight in occurrence estimates. This typically results in overestimation of occurrence rates.

This observer-excitement bias is very difficult to correct for, especially when combining blind surveys with different observing strategies.

Here, we present a long-term observational effort to study the OGP-ILP (anti-)correlation by probing the ILP occurrence in systems with and without OGP. We constructed two comparative samples of stars, with and without OGP, ensuring the homogeneity of the two samples in terms of stellar properties and observing strategy, to avoid any bias that could affect the comparison. In this first article of a series, we focus on the sample definitions and on the first three detections (two new and one confirmation) of ILP in OGP hosts around HD 23079, HD 196067, and HD 86226, based on the first two years of observations of our program. These three targets were known to host an OGP (see Tinney et al. 2002 for HD 23079, Marmier et al. 2013 for HD 196067, and Arriagada et al. 2010 for HD 86226). Moreover, a transiting ILP has been detected by TESS (Transiting Exoplanet Survey Satellite) and confirmed with follow-up PFS (Planet Finder Spectrograph) RV measurements (Teske et al. 2020) for HD 86226.

In Sect. 2, we present the two samples and the observing strategy. In Sect. 3, we present the properties of the three hosts, and in Sect. 4 we describe the observations. We present our analysis methods in Sect. 5, and our results for the three targets in Sect. 6. Finally, in Sect. 7, we present preliminary conclusions from these detections. A statistical analysis of the full sample to derive the underlying occurrences of ILP with and without OGP will be performed in subsequent articles.

## 2. Sample definitions and observing strategy

We constructed our samples from two historical volume-limited giant-planet search surveys conducted with CORALIE (Udry et al. 2000) and HARPS (Lo Curto et al. 2010) spectrographs. These two historical surveys comprise a total of 2487 main-sequence G and K dwarfs, volume-limited to 57 pc, and were conducted at a typical precision of about 5 m/s for CORALIE and 2–3 m/s for HARPS over more than two decades.

We selected less active stars ($\log R'_{HK} < -4.6$) and considered only stars closer than 52 pc. We defined two subsamples for our comparative study, composed of stars with similar mass distributions ranging from 0.65 to $1.4 m_⊕$:

1. “Giant” sample (all 57 systems matching the criteria):

- with at least one giant planet in their outer region ($m \sin i > 100 m_⊕$, $a > 0.5$ AU),

- no giant planet in their inner region ($a < 0.5$ AU),

- no highly eccentric orbits ($e < 0.7$),

- no high-mass companions ($m < 30 m_⊕$ at 5 AU) based on Hipparcos-Gaia proper motion anomalies (PMA, see Kervella et al. 2022),

2. “Flat” sample (55 systems randomly selected among 73 systems matching the criteria):

- at least 25 HARPS and 25 CORALIE pre-existing measurements spread over at least 8 years,
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Masses (top panel) and metallicities (bottom panel) of the stars in our Giant and Flat samples (see Sect. 2).</figcaption></figure>

- no signature of a giant planet from the analysis of available RVs,

- no signature of a giant planet from the analysis of Hipparcos-Gaia PMA.

We obtained two successive ESO Large Programmes (LPs) on HARPS (108.22KV) and HARPS+ESPRESSO (112.25YG). In both campaigns, the aim was to explore the inner parts of these systems by taking advantage of preexisting measurements. We searched for ILPs in the super-Earth (SE; 3–10 $m_⊕$) and Neptune-like (NE; 10–30 $m_⊕$) ranges, with periods shorter than 50 d.

The goal of the first LP (108.22KV) was to obtain approximately 75 HARPS measurements for 26 systems belonging to the Giant sample. The second LP (112.25YG), which is on-going at the time of writing, aims to gather 50 HARPS measurements for all 112 systems and 25 ESPRESSO measurements for half of the sample (25 systems from the Giant sample and 27 from the Flat sample). Most systems already had HARPS measurements, with some having reached the 50, or even 75, measurement thresholds, allowing reliance on archival data where appropriate. We used exposure times of 15 minutes for HARPS and 12 minutes for ESPRESSO (split in sub-exposures for bright targets) to average out stellar oscillations (Dumusque et al. 2011).

We note that the two historical giant-planet search surveys were subject to the observer-excitement bias described above. Before our search for inner planets began, the stars in the Giant sample had a median of 109 CORALIE measurements, while the stars in the Flat sample had a median of 34 CORALIE measurements. However, the detection capacity for low-mass inner planets in our samples is dominated by the high-precision HARPS and ESPRESSO measurements, for which we ensured a homogeneous observing strategy. Moreover, these historical

surveys are blind-search surveys (i.e., not transit follow-ups), so there is no bias toward stars with known inner companions. Thus, the observer-excitement bias in giant-planet search surveys could have increased the probability of detecting additional OGPs in systems that show a hint of OGP presence. However, the effect of this bias on the OGP–ILP correlation should be limited.

Fig. 1 shows the distribution of masses and metallicities of our two samples, as well as those of the sub-samples observed with ESPRESSO. These masses and metallicities were computed using the ARES+MOOG method (Sousa et al. 2021, see also Sect. 3). The mass distributions of the two samples are very similar (Fig. 1). Regarding metallicities, the Giant sample is biased toward higher metallicities than the Flat sample, due to the metallicity-planet correlation (Gonzalez 1998; Santos et al. 2001). However, both samples contain a large number of targets in the range −0.4 < [Fe/H] < 0.2, allowing an unbiased comparison (with respect to metallicity) by excluding the most extreme targets in each sample. In the longer term, the Flat sample may be slightly refined for future observational programs to exclude the most metal-poor stars and include additional metal-rich stars that were randomly excluded.

Here, we present three inner planet detections from our analysis of the 26 systems observed during the first LP (108.22KV). Global statistical analysis of the sample will be performed in subsequent articles of this series.

### 3. Stellar properties

This section describes the fundamental parameters of HD 23079, HD 196067, and HD 86226, which are part of our Giant sample. Obtaining accurate stellar masses is crucial, as it directly impacts our estimates of the planetary (minimum) masses inferred from the RV time series. Values of effective temperature $T_{eff}$, stellar mass $M_*$, stellar luminosity $L_*$, stellar radius $R_*$, and distance $d$ were estimated using the ARIADNE (spectral energy distribution bAyesian model) fitting Er) Python package (Vines & Jenkins 2022).

The ARIADNE package fits the spectral energy distribution (SED) of a star using multiple stellar atmosphere model grids, including PHOENIX V2 (Husser et al. 2013), BT-Settl, BT-NextGen, BT-Cond (Allard et al. 2012), Castelli & Kurucz (Castelli & Kurucz 2003), and Kurucz (Kurucz 1993). The package automatically retrieves available broadband photometric data from various sources and performs an SED fit using the different models. The resulting fits are then compared using Bayesian model averaging, which yields estimates of the stellar parameters. All of the results presented here were obtained by running ARIADNE with default priors (see Table 2 of Vines & Jenkins 2022).

Because the SED does not precisely constrain the stellar metallicity [Fe/H] and surface gravity $g$, we derived these values using the ARES+MOOG method (Sneden 1973; Sousa et al. 2007, 2008, 2015, 2021; Sousa 2014). Additionally, surface gravity values were obtained from the Gaia Early Data Release 3 (Gaia Collaboration 2016; Riello et al. 2021) parallaxes and photometry. We computed the log $R'_{HK}$ time series using Gomes da Silva et al. (2021); Rutten (1984) and derived estimates of the rotation period and age using Mamajek & Hillenbrand (2008). A summary of the derived stellar parameters is given in Table 1.

HD 23079 is an F9.5 V dwarf at a distance $d = 33.5$ pc with a magnitude of $V = 7.11$. The ARIADNE fit yields an effective temperature $T_{eff} = 5994$ K and a mass estimate
<table><caption>**Table 1.** Stellar parameters for HD 23079, HD 196067, and HD 86226.</caption><thead><tr><th>Star</th><th>HD 23079</th><th>HD 196067</th><th>HD 86226</th><th></th></tr></thead><tbody><tr><td>Spectral Type</td><td>F9.5V</td><td>G0V</td><td>G2V</td><td>a, b</td></tr><tr><td>π (mas)</td><td>29.8633 ± 0.0194</td><td>25.0327 ± 0.0205</td><td>21.9301 ± 0.0267</td><td>c</td></tr><tr><td>V (mag)</td><td>7.11</td><td>6.44</td><td>7.93</td><td>a</td></tr><tr><td>B (mag)</td><td>7.68</td><td>7.05</td><td>8.563</td><td>a</td></tr><tr><td>d (pc)</td><td>33.50<sup>+0.05</sup><sub>-0.02</sub></td><td>40.0<sup>+0.05</sup><sub>-0.04</sub></td><td>45.68<sup>+0.09</sup><sub>-0.08</sub></td><td>d</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5994<sup>+55</sup><sub>-64</sub></td><td>6072<sup>+83</sup><sub>-76</sub></td><td>6007<sup>+75</sup><sub>-71</sub></td><td>d</td></tr><tr><td>[Fe/H] (dex)</td><td>-0.129 ± 0.041</td><td>0.213 ± 0.042</td><td>-0.015 ± 0.041</td><td>e</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.37 ± 0.02</td><td>4.09 ± 0.02</td><td>4.41 ± 0.02</td><td>e, c</td></tr><tr><td>R<sub>*</sub>(R<sub>☉</sub>)</td><td>1.08<sup>+0.02</sup><sub>-0.03</sub></td><td>1.71<sup>+0.03</sup><sub>-0.04</sub></td><td>1.03<sup>+0.02</sup><sub>-0.02</sub></td><td>d</td></tr><tr><td>M<sub>*</sub>(M<sub>☉</sub>)</td><td>1.16<sup>+0.02</sup><sub>-0.03</sub></td><td>1.16<sup>+0.13</sup><sub>-0.09</sub></td><td>1.08<sup>+0.05</sup><sub>-0.05</sub></td><td>d</td></tr><tr><td>L<sub>*</sub>(L<sub>☉</sub>)</td><td>1.35<sup>+0.08</sup><sub>-0.08</sub></td><td>3.54<sup>+0.27</sup><sub>-0.20</sub></td><td>1.24<sup>+0.09</sup><sub>-0.07</sub></td><td>d</td></tr><tr><td>log R'<sub>HK</sub></td><td>-4.929 ± 0.014</td><td>-5.033 ± 0.022</td><td>-4.946 ± 0.012</td><td>f</td></tr><tr><td>P<sub>rot.</sub> (d)</td><td>15.0 ± 2.6</td><td>22.5 ± 2.9</td><td>22.8 ± 3.0</td><td>g</td></tr><tr><td>Age (Gyr)</td><td>2.6 ± 1.2</td><td>4.2 ± 2.0</td><td>3.8 ± 1.8</td><td>g</td></tr></tbody></table>

**Notes.** The uncertainties in this table only reflect the formal error propagation. For some parameters (e.g., the masses), the systematic errors in stellar models dominate these formal uncertainties (e.g., Tayar et al. 2022). (a) Gray et al. (2006); (b) Houk & Smith-Moore (1988); (c) Riello et al. (2021); (d) Obtained from ARIADNE (Vines & Jenkins 2022); (e) Obtained using the ARES+MOOG method (Sousa et al. 2021); (f) Gomes da Silva et al. (2021); Rutten (1984); (g) Mamajek & Hillenbrand (2008).

$M_* = 1.16 M_☉$. The star has a surface gravity $log g = 4.37$ cm s$^{-2}$ and a metallicity $[Fe/H] = -0.129$ dex.

HD 196067 is a G0V star with magnitude $V = 6.44$. From the ARIADNE analysis, we find that this star is located at a distance $d = 40$ pc, has an effective temperature $T_{eff} = 6072$ K, and a mass $M_* = 1.16 M_☉$. The surface gravity and metallicity are $log g = 4.09$ cm s$^{-2}$ and $[Fe/H] = 0.213$ dex, respectively.

HD 86226 is a G2V dwarf with magnitude $V = 7.93$ located at a distance $d = 45.68$ pc. From the SED fit, we obtain an effective temperature $T_{eff} = 6007$ K. This star has an estimated mass $M_* = 1.08 M_☉$, a metallicity $[Fe/H] = -0.015$ dex, and a surface gravity $log g = 4.41$ cm s$^{-2}$.

## 4. Observations

This section describes the combined RV time series for the three targets from CORALIE and HARPS observations, with the addition of PFS observations for HD 86226.

The CORALIE instrument is an echelle spectrograph installed on the Euler Swiss telescope at La Silla Observatory in Chile since 1998. It underwent upgrades in 2007 and 2014, which led to possible RV offsets. We consider the three versions of CORALIE as different instruments: CORALIE98, CORALIE07, and CORALIE14. Accordingly, we adjusted for a different offset and a different jitter term for each CORALIE version. The HARPS instrument is also an echelle spectrograph at La Silla Observatory, installed on the 3.6m telescope since 2003. We distinguish two versions of HARPS as different instruments, HARPS03 and HARPS15, corresponding to before and after a fiber link upgrade performed in May 2015 (Lo Curto et al. 2015).

For HD 86226, we also used the PFS RV measurements from Teske et al. (2020). The PFS instrument is an echelle spectrograph installed on the 6.5 m Magellan II telescope at Las Campanas Observatory in Chile. Its CCD was upgraded in 2018, and we distinguish two versions of the instrument, PFS1 and PFS2, corresponding to before and after this upgrade. We used these PFS RVs in addition to the CORALIE and HARPS time series.

<table><caption>**Table 2.** Summary of observations used in our analysis.</caption><thead><tr><th></th><th></th><th colspan="2">HD 23079</th><th colspan="2">HD 196067</th><th colspan="2">HD 86226</th></tr><tr><th></th><th></th><th>#</th><th>Avg. err.</th><th>#</th><th>Avg. err.</th><th>#</th><th>Avg. err.</th></tr></thead><tbody><tr><td rowspan="2">CORALIE</td><td>Spectra</td><td>129</td><td>3.89</td><td>116</td><td>4.02</td><td>116</td><td>4.41</td></tr><tr><td>Nights</td><td>123</td><td>3.88</td><td>112</td><td>3.83</td><td>114</td><td>4.34</td></tr><tr><td rowspan="2">HARPS</td><td>Spectra</td><td>77</td><td>0.67</td><td>76</td><td>0.66</td><td>74</td><td>0.84</td></tr><tr><td>Nights</td><td>76</td><td>0.67</td><td>75</td><td>0.63</td><td>73</td><td>0.83</td></tr><tr><td rowspan="2">PFS</td><td>Spectra</td><td>-</td><td>-</td><td>-</td><td>-</td><td>105</td><td>1.4</td></tr><tr><td>Nights</td><td>-</td><td>-</td><td>-</td><td>-</td><td>59</td><td>1.34</td></tr></tbody></table>

A summary of the observations used in the analysis of the three targets is given in Table 2. The table lists both the number of spectra and the number of unique nights of observations for each instrument. In the following, we used nightly binned time series.

## 5. Analysis

We performed two analyses for the three systems: a simple periodogram analysis (Sect. 5.1) and posterior inference using a nested sampling approach (Sect. 5.2).

### 5.1. Periodogram analysis

As an initial, simple, and fast analysis, we performed a period search using the open-source kepmodel<fnref n="2" /> Python package (Delisle et al. 2016, 2020). We first defined a base model of the RV time series, with an offset and a jitter term for each

<footnote n="2">² Available at: https://gitlab.unige.ch/delisle/kepmodel</footnote>
instrument. The jitter term is added quadratically to the RV uncertainties. We performed a maximum-likelihood fit of the model, restricting the jitter terms to be lower than 20 m s⁻¹ (see Table 3). We then searched sequentially for planetary signals in the residuals of the RV time series by computing a periodogram and the associated false alarm probability (FAP; see Baluev 2008; Delisle et al. 2020). For low FAP values (below 0.1%), we selected the highest peak in the periodogram and computed an estimate of the planet's orbital parameters (see Delisle et al. 2016; Delisle & Ségransan 2022). We then performed a global maximum likelihood fit, adjusting all planetary parameters as well as the instruments' offsets and jitter terms. We iterated the procedure until no significant periodic signal was found in the residuals (i.e., FAP > 0.1%).

## 5.2. Posterior sampling

To better determine the orbital architecture of the systems and estimate uncertainties on the orbital parameters, we explored the parameter space using a nested sampling approach. The same model was used for the three targets, but with different prior distributions for the outer companions. Since the detection of these outer companions is not in question, we set some priors based on the previous periodogram analysis. This choice does not significantly impact the resulting posterior distributions. For the inner part of the systems, we adopted weakly informative priors, since the goal was to detect lower-mass planets, if present.

The model then consists of the following components: a Keplerian function for the outer companion, up to two additional Keplerians for the inner part of the system, between-instrument RV offsets, a constant systemic velocity, and a Gaussian noise model per instrument that includes individual jitter terms. The prior distributions for the parameters are listed in Table 3.

We used the Gaia DR3 radial velocity measurement (Katz et al. 2023) to set the prior on the systemic velocity, $v_{sys}$, and the associated uncertainty to assign a Gaussian prior for the between-instrument RV offsets. Even if relatively informative, these priors are much wider than the posteriors and do not influence the results. For the jitter terms, we assigned a modified log-uniform prior (see e.g. Gregory 2005) with a knee at 1 m s⁻¹ and an upper limit of 20 m s⁻¹. The number of additional Keplerians, $N_p$, was free in this analysis and was assigned a uniform prior from 0 to 2, motivated by the fact that the periodogram analysis from Section 5.1 showed only one inner companion in each system. The orbital parameters of these two Keplerians share the same priors: a log-uniform distribution for the periods, between 1 day and the time span of the time series, a modified log-uniform distribution for the semi-amplitudes (with a knee at 1 m s⁻¹ and an upper limit of 100 m s⁻¹), a Kumaraswamy distribution (Kumaraswamy 1980) for the eccentricities, and uniform priors for the argument of periastron, $\omega$, and the mean anomaly, $M_0$. The Kumaraswamy distribution is very similar to the beta distribution proposed by Kipping (2013), but is easier to implement in nested sampling algorithms since its inverse cumulative distribution function has a closed form (see also Stevenson et al. 2025).

We sampled from the posterior distribution using kima (Faria et al. 2018), a package<fnref n="3" /> dedicated to the Bayesian analysis of RV time series. The code uses the diffusive nested sampling algorithm (Brewer et al. 2011) to explore the posterior and calculate the marginal likelihood, or evidence, of the model. The algorithm is a variant of classic nested sampling

<table><caption>**Table 3.** Prior distributions used in the RV analysis.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Prior</th></tr></thead><tbody><tr><td colspan="3"><strong>General</strong></td></tr><tr><td><i>v</i><sub>sys</sub></td><td>m s<sup>-1</sup></td><td><i>G</i> (RV<sub>Gaia</sub>, σRV<sub>Gaia</sub>)</td></tr><tr><td><i>v</i><sub>offset,i</sub></td><td>m s<sup>-1</sup></td><td><i>G</i> (0, σRV<sub>Gaia</sub>)</td></tr><tr><td><i>j</i></td><td>m s<sup>-1</sup></td><td><i>MLU</i> (1, 20)</td></tr><tr><td colspan="3"><strong>Outer companion</strong></td></tr><tr><td><i>P</i></td><td>days</td><td><i>LU</i> (0.8 &middot; <i>P</i><sub>guess</sub>, 1.2 &middot; <i>P</i><sub>guess</sub>)</td></tr><tr><td><i>K</i></td><td>m s<sup>-1</sup></td><td><i>U</i> (0.5 &middot; <i>K</i><sub>guess</sub>, 1.5 &middot; <i>K</i><sub>guess</sub>)</td></tr><tr><td><i>e</i></td><td></td><td><i>U</i> (0, 1)</td></tr><tr><td><i>ω</i></td><td>rad</td><td><i>U</i> (0, 2π)</td></tr><tr><td><i>M</i><sub>0</sub></td><td>rad</td><td><i>U</i> (0, 2π)</td></tr><tr><td colspan="3"><strong>Additional Keplerians</strong></td></tr><tr><td><i>N</i><sub>p</sub></td><td>days</td><td><i>U</i><sub>i</sub> (0, 2)</td></tr><tr><td><i>P</i></td><td>days</td><td><i>LU</i> (1, Δ<i>t</i>)</td></tr><tr><td><i>K</i></td><td>m s<sup>-1</sup></td><td><i>MLU</i> (1, 100)</td></tr><tr><td><i>e</i></td><td></td><td><i>K</i> (0.867, 3.03)</td></tr><tr><td><i>ω</i></td><td>rad</td><td><i>U</i> (0, 2π)</td></tr><tr><td><i>M</i><sub>0</sub></td><td>rad</td><td><i>U</i> (0, 2π)</td></tr></tbody></table>

**Notes.** The distributions are abbreviated as follows: *G* = Gaussian with mean and standard deviation; *MLU* = modified log-uniform with knee and upper limit; *LU* = log-uniform with lower and upper limits; *U* = uniform with lower and upper limits; *K* = Kumaraswamy with parameters *a* and *b*. *P<sub>guess</sub>* and *K<sub>guess</sub>* denote the maximum likelihood estimates for the period and semi-amplitude of the outer companion. Δ<i>t</i> is the time span of the time series.

(Skilling 2006) with some unique advantages. It explores a mixture of nested distributions, each successively occupying smaller regions within the prior, guided by the likelihood function. It can handle distributions showing multimodality, strong parameter correlations, and phase transitions (see e.g. Buchner 2023), much like those typically encountered in the exoplanet detection problem (e.g. Ford 2005; Gregory 2005; Feroz et al. 2011). Moreover, the algorithm can handle so-called trans-dimensional models (Brewer 2014; Brewer & Donovan 2015), in which the number of parameters can change depending on an unknown number of components (Keplerian functions, in our case). We ran the algorithm for a fixed number of steps (500 000), which guaranteed a sufficient effective sample size (ESS).

# 6. Results

Fig. 2 shows successive periodograms of the raw RV time series and after subtracting planetary signal for the three targets (HD 23079, HD 196067, and HD 86226). In all three systems, the periodogram analysis clearly indicates the detection of two planets, an OGP and an ILP. These detections are significant, as the FAP of all six planet signals are lower than 10⁻⁷. The FAP of the residuals after subtracting the planetary signals are above 10% in all three time series, indicating that no clear signal remains unmodeled in the residuals.

We show the true inclusion probability (TIP) periodograms based on the kima posterior sampling in Fig. 3. The TIP (see Hara et al. 2022) is the posterior probability (conditioned on the data) of having at least one planet in a small period interval around a given period *P*, or equivalently, in a small frequency interval [*f* − Δ*f*/2, *f* + Δ*f*/2] around a frequency *f* = 1/*P*. The TIP has been mathematically shown to be an optimal exoplanet

<footnote n="3">³ Available at: https://github.com/kima-org/kima</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Periodograms of the raw RV time series (top panel) as well as of the residuals after subtracting the OGP (middle panel) and the ILP (bottom panel), for HD 23079 (**left**), HD 196067 (**center**), and HD 86226 (**right**).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Posterior distributions for the orbital periods from the analysis of the three targets (**left**: HD 23079, **center**: HD 196067, **right**: HD 86226). The plots show the true inclusion probability (TIP), highlighting the peaks corresponding to the detected inner planets. Because the Keplerian signal of the outer planet is always included in the model, the TIP value for the outer planet period is, by definition, unity.</figcaption></figure>

detection criterion, in the sense that it maximizes the expected number of true detections for a given tolerance to false detections (see Hara et al. 2024). The TIP periodograms of Fig. 3 were obtained using $\Delta f = 1/T_{\text{span}}$ (where $T_{\text{span}}$ is the time span of the observations) and by scanning the frequency *f* over a finer grid ($\Delta f/10$) to ensure that local maxima of the TIP were not missed. The TIP is marginalized over the number of planets included in the model. In the *kima* runs, we always included the OGP in the model and allowed the number of additional planets to vary between zero and two (thus the total number of planets is free between one and three). Therefore, by construction, the TIP is one at the OGP's period because we enforced the presence of this already known planet in the model. For the ILP, the TIP values exceed 99% for all three targets (see Fig. 3), which confirms the detection of these planets. As in the classical periodograms, we do not find any additional significant peak in the TIP periodograms.

The parameters of the maximum-likelihood solution and of the posterior sampling are given in Table A.1 for all three systems. Table 4 provides a shorter version of this table that contains only the posteriors of the orbital parameters. The two estimates are in good agreement for all parameters. As final estimates,

we quote the maximum-likelihood values and the 68% credible intervals from the posteriors. The RV time series and the residuals from the maximum-likelihood solution are shown in Fig. 4 and are phase-folded in Fig. 5. Our solution for HD 86226 agrees with Teske et al. (2020), with a slight refinement to the mass estimate for the inner planet (c).

From the *kima* analysis, we also derived so-called compatibility limits, corresponding to the orbital periods and masses of undetected planets that remain compatible with the available data (see, e.g., Standing et al. 2022 and Figueira et al. 2025). In practice, these limits were obtained from all posterior samples with $N_p = 2$ (thus with a total of three Keplerians in the model, accounting for the OGP) and are therefore available from a single run of the algorithm. These limits do not rely on injection and recovery of signals (sinusoidal or otherwise), do not require the subtraction of a best-fit solution, and marginalize over other orbital parameters, such as eccentricity. The compatibility limits for the three targets are shown in Fig. 6. The adjusted HARPS jitter is larger for HD 196067 than for the other two targets (see Table A.1), due to a larger scatter in the residuals for this star. This results in slightly worse detection limits for this source (see Fig. 6). For HD 86226, the compatibility limits shown in
<table><caption>**Table 4.** Posteriors of the orbital parameters for HD 23079, HD 196067, and HD 86226.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>HD 23079</th><th>HD 196067</th><th>HD 86226</th></tr></thead><tbody><tr><td>P<sub>b</sub></td><td>d</td><td>735.7<sup>+0.4</sup><sub>-0.4</sub></td><td>3397<sup>+17</sup><sub>-14</sub></td><td>1616<sup>+46</sup><sub>-28</sub></td></tr><tr><td>U<sub>0,b</sub> = M<sub>0,b</sub> + ω<sub>b</sub><sup>a</sup></td><td>deg</td><td>107.0<sup>+0.9</sup><sub>-0.9</sub></td><td>-53.0<sup>+2.6</sup><sub>-2.2</sub></td><td>-66<sup>+16</sup><sub>-12</sub></td></tr><tr><td>K<sub>b</sub></td><td>m s<sup>-1</sup></td><td>54.4<sup>+0.3</sup><sub>-0.4</sub></td><td>80<sup>+3</sup><sub>-4</sub></td><td>7.8<sup>+0.8</sup><sub>-0.6</sub></td></tr><tr><td>e<sub>b</sub></td><td></td><td>0.072<sup>+0.012</sup><sub>-0.011</sub></td><td>0.561<sup>+0.011</sup><sub>-0.012</sub></td><td>0.07<sup>+0.09</sup><sub>-0.05</sub></td></tr><tr><td>ω<sub>b</sub><sup>b</sup></td><td>deg</td><td>29<sup>+15</sup><sub>-7</sub></td><td>142.0<sup>+1.5</sup><sub>-1.8</sub></td><td>37<sup>+100</sup><sub>-180</sub></td></tr><tr><td>m<sub>b</sub> sin i<sub>b</sub><sup>c</sup></td><td>M<sub>J</sub></td><td>2.66<sup>+0.05</sup><sub>-0.05</sub></td><td>5.42<sup>+0.20</sup><sub>-0.26</sub></td><td>0.49<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>P<sub>c</sub></td><td>d</td><td>5.74894<sup>+0.00022</sup><sub>-0.00021</sub></td><td>4.60111<sup>+0.00024</sup><sub>-0.00024</sub></td><td>3.98450<sup>+0.00023</sup><sub>-0.00022</sub></td></tr><tr><td>U<sub>0,c</sub> = M<sub>0,c</sub> + ω<sub>c</sub><sup>a</sup></td><td>deg</td><td>156<sup>+10</sup><sub>-10</sub></td><td>-1<sup>+18</sup><sub>-17</sub></td><td>-10<sup>+24</sup><sub>-22</sub></td></tr><tr><td>K<sub>c</sub></td><td>m s<sup>-1</sup></td><td>2.65<sup>+0.26</sup><sub>-0.27</sub></td><td>3.6<sup>+0.4</sup><sub>-0.4</sub></td><td>2.82<sup>+0.25</sup><sub>-0.26</sub></td></tr><tr><td>e<sub>c</sub></td><td></td><td>0.07<sup>+0.09</sup><sub>-0.05</sub></td><td>0.12<sup>+0.12</sup><sub>-0.09</sub></td><td>0.05<sup>+0.07</sup><sub>-0.04</sub></td></tr><tr><td>ω<sub>c</sub><sup>b</sup></td><td>deg</td><td>247<sup>+64</sup><sub>-166</sub></td><td>278<sup>+41</sup><sub>-104</sub></td><td>170<sup>+108</sup><sub>-94</sub></td></tr><tr><td>m<sub>c</sub> sin i<sub>c</sub><sup>c</sup></td><td>M<sub>⊕</sub></td><td>8.1<sup>+0.8</sup><sub>-0.8</sub></td><td>10.1<sup>+1.2</sup><sub>-1.2</sub></td><td>7.7<sup>+0.7</sup><sub>-0.7</sub></td></tr><tr><td>Epoch</td><td>BJD</td><td>2 455 000</td><td>2 455 000</td><td>2 455 000</td></tr></tbody></table>

**Notes.** An extended version of this table with the maximum likelihood estimates and kima posteriors of all model parameters is provided in Table A.1. (a) Mean argument of latitude at the reference epoch ($U_0 = M_0 + \omega = \lambda_0 - \Omega$, where $M_0$, $\lambda_0$, $\omega$, and $\Omega$ are the mean anomaly, mean longitude, argument of periastron, and longitude of the ascending node, respectively). (b) Argument of periastron of the stellar orbit. (c) The uncertainties on the planetary masses take into account the stellar mass formal uncertainties, but do not include systemic errors from stellar models (see Table 1).

Fig. 6 include PFS RVs in addition to CORALIE and HARPS data. However, excluding PFS data does not significantly change the detection limits for orbital periods shorter than a year (see Appendix B).

All three targets were observed by TESS. We checked whether the inner planets detected for each star transit in the available TESS data. This analysis is presented in Appendix C. We conclude that HD 23079 and HD 196067 are not transiting, while the known transit of HD 86226 is clearly recovered (see Teske et al. 2020).

To constrain the inclinations of the OGP and thereby their true masses, we jointly fit RV and astrometric data from the HIPPARCOS-Gaia Catalog of Accelerations (HGCA; Brandt 2021) using orvara (Brandt et al. 2021). The HGCA provides the $\chi^2$ value for each star, which should follow a $\chi^2$ distribution with 2 degrees of freedom if the star lacks acceleration. The $\chi^2$ is 0.709 (p-value of 70%) for HD 23079, 64.9 (p-value of $8 \times 10^{-15}$) for HD 196067, and 1.12 (p-value of 57%) for HD 86226. Thus, we detect a highly significant acceleration for HD 196067, and nonsignificant accelerations for HD 23079 and HD 86226. We used a customized version of orvara$^4$, which accommodates a larger set of priors, in particular for between-instrument RV offsets. All priors were set as in Sect. 5.2, except for RV jitters (log-uniform distribution from 0.1 m s$^{-1}$ to 20 m s$^{-1}$), and for the semi-major axis (log-uniform from 0.0025 AU to 500 AU) and mass (log-uniform from $10^{-7} M_\odot$ to $1 M_\odot$) of the planets, which were explored instead of the period and RV semi-amplitude. We included a Gaussian prior on the period of the ILP (the maximum likelihood solution from Table A.1), to facilitate the convergence

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Time series (top panel) and residuals (bottom panel) of the maximum-likelihood fit for each target (HD 23079, HD 196067, and HD 86226).</figcaption></figure>

of the fit. We excluded the binary companion to HD 196067, HD 196068, from the orvara fit, because it lies at a separation of >16 arcseconds (≈ 800 AU, Gaia Collaboration 2023), too far to be detected in our RVs or in HGCA astrometry of HD 196067.

We ran orvara with ten temperatures, each with 150 walkers, for 150 000 steps. As expected, we are unable to constrain the inclination of the three ILPs. We also find no inclination constraints for the OGP of HD 23079 and HD 86226, which is also

<footnote>$^4$ Available at: https://gitlab.unige.ch/Damien.Segransan/orvara</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Top panel: residuals of the maximum-likelihood fit phase-folded at the periods of the two detected planets, showing the outer giant (left) and lower-mass inner planet (right) for each target (HD 23079, HD 196067, and HD 86226). Bottom panel: same residuals, but showing only the HARPS data.</figcaption></figure>

expected given the nonsignificant Hipparcos-Gaia accelerations for these two stars.

The posteriors of the *orvara* orbital fit for HD 196067 are shown in Fig. 7. We find that the true mass of the OGP of HD 196067 is $m_b = 11.4 \pm 1.6 M_J$. The inclination distribution is bimodal, with two almost symmetric modes at $i = 29.7_{-3.8}^{+5.0}$ deg and $i = 149.2_{-5.1}^{+3.9}$ deg. The mode at 29.7 deg appears slightly favored by the data. We find a semi-major axis of $4.71_{-0.10}^{+0.11}$ AU.

These results update those presented in Li et al. (2021), providing finer constraints on the orbital parameters of the OGP (a factor of two on the semi-major axis) thanks to our additional RV data.

## 7. Discussion

This article presents the three inner planet detections obtained by monitoring a sample of 26 stars with known outer giant planets
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Compatibility limits from the analysis of the three targets (**left**: HD 23079, **center**: HD 196067, **right**: HD 86226). The plots show a color map of the posterior density for planet mass in orbital period bins, for all samples with $N_p = 2$ (i.e., with a total of three Keplerians, including the OGP). Green points indicate detected inner planets with corresponding mass uncertainties. Dashed and dotted blue lines indicate masses corresponding to 1 and 5 m/s semi-amplitude Keplerians, respectively. Masses are calculated using the stellar masses from Table 1. Red boxes highlight the SE (3–10 $m_⊕$) and NE (10–30 $m_⊕$) mass ranges within 50 d.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Posterior distributions from the *orvara* fit of HD 196067 (see Sect. 6).</figcaption></figure>

using HARPS. A dedicated statistical study including these 26 stars and additional well-monitored stars with and without outer giant planets will be presented in a subsequent paper. We note here that two of these inner planets lie in the SE range (3–10 $m_⊕$), and one is in the NE range (10–30 $m_⊕$), slightly above the SE-NE limit.

The estimated occurrences of SE and NE below 50 d in the general population are 0.11 ± 0.03 and 0.08 ± 0.02, respectively (Bashi et al. 2020). If the occurrence rates of SE and NE were the same in our sample as in the general population (i.e., without impact from the giant planets), we would typically expect to detect 2.9 SE and 2.1 NE (about five inner planets in total). From the compatibility limits of the three targets presented here (see Fig. 6), we see that our observing strategy enables the detection

of NE below 50 d with very high completeness. For SE, completeness is very poor above 10 d, and we aim to improve this in the long term with additional HARPS and ESPRESSO measurements (on-going ESO LP). Since we used the same observing strategy for the 26 targets in our sample, we can extrapolate, as a first approximation, that similar detection limits are reached for the 23 other targets, for which we do not find any ILP.

We detect fewer inner planets than expected from the general population occurrence rates, but given the small-number statistics and without a proper statistical analysis correcting for completeness, we cannot conclude on a slight enhancement or inhibition of ILP in the presence of OGP. However, these results are difficult to reconcile with the claim of Zhu & Wu (2018); Bryan et al. (2019) that $p(\text{SE} \cup \text{NE}|\text{OGP})$ is close to 90%.

Interestingly, the three ILPs detected in this study all have periods below 6 d. While this could be due to small-number statistics and reduced detection efficiency at longer periods (see Fig. 6), these three detections are significant, and we should have been able to detect, albeit with lower confidence, similar-mass planets up to 20 d. This observation, if confirmed with a rigorous statistical analysis, would be in agreement with the findings of Van Zandt et al. (2025) that ILPs tend to be found at shorter periods in systems with an OGP than in the general population.

## Data availability

The RV data for the three targets discussed in this work, as well as *kima* posterior samples, are available on the DACE platform: https://doi.org/10.82180/dace-inner001.

*Acknowledgements.* Based on observations collected at the European Southern Observatory under ESO programmes 072.C-0488, 192.C-0852, 0101.C-0275, 0103.C-0785, 1102.C-0923, and 108.22KV. The authors acknowledge the financial support of the Swiss National Science Foundation (SNSF), supported since May 2022 over the grant 200020_205010. The 120 cm EULER telescope and the CORALIE spectrograph were funded by the SNSF and the University of Geneva. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the SNSF under grants 51NF40_182901 and 51NF40_205606. 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 dace.unige.ch. This
research has made use of the SIMBAD database, CDS, Strasbourg Astronomical Observatory, France (Wenger et al. 2000). XD 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). NCS is Funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020 DOI: 10.54499/UIDB/04434/2020, UIDP/04434/2020 DOI: 10.54499/UIDP/04434/2020.

## References

Allard, F., Homeier, D., & Freytag, B. 2012, Philos. Trans. Roy. Soc. Lond. A, 370, 2765

Arriagada, P., Butler, R. P., Minniti, D., et al. 2010, ApJ, 711, 1229

Baluev, R. V. 2008, MNRAS, 385, 1279

Barbato, D., Sozzetti, A., Desidera, S., et al. 2018, A&A, 615, A175

Bashi, D., Zucker, S., Adibekyan, V., et al. 2020, A&A, 643, A106

Bonomo, A. S., Dumusque, X., Massa, A., et al. 2023, A&A, 677, A33

Bonomo, A. S., Naponiello, L., Pezzetta, E., et al. 2025, A&A, 700, A126

Brandt, T. D. 2021, ApJS, 254, 42

Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021, AJ, 162, 186

Brewer, B. J. 2014, arXiv e-prints [arXiv:1411.3921]

Brewer, B. J., & Donovan, C. P. 2015, MNRAS, 448, 3206

Brewer, B. J., Pártay, L. B., & Csányi, G. b. 2011, Stat. Comput., 21, 649

Bryan, M. L., & Lee, E. J. 2024, ApJ, 968, L25

Bryan, M. L., & Lee, E. J. 2025, ApJ, 982, L7

Bryan, M. L., Knutson, H. A., Lee, E. J., et al. 2019, AJ, 157, 52

Buchner, J. 2023, Stat. Surv., 17, 169

Castelli, F., & Kurucz, R. L. 2003, Proceedings of the IAU Symp., 210, Modelling of Stellar Atmospheres, eds. N. Piskunov et al., A20

Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444

Delisle, J. B., & Ségransan, D. 2022, A&A, 667, A172

Delisle, J.-B., Ségransan, D., Buchschacher, N., & Alesina, F. 2016, A&A, 590, A134

Delisle, J. B., Hara, N., & Ségransan, D. 2020, A&A, 635, A83

Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140

Faria, J. P., Santos, N. C., Figueira, P., & Brewer, B. J. 2018, JOSS, 3, 487

Feroz, F., Balan, S. T., & Hobson, M. P. 2011, MNRAS, 415, 3462

Figueira, P., Faria, J. P., Silva, A. M., et al. 2025, A&A, 700, A174

Ford, E. B. 2005, AJ, 129, 1706

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Gomes da Silva, J., Santos, N. C., Adibekyan, V., et al. 2021, A&A, 646, A77

Gonzalez, G. 1998, A&A, 334, 221

Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161

Gregory, P. C. 2005, ApJ, 631, 1198

Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022, A&A, 663, A14

Hara, N. C., de Poyferré, T., Delisle, J.-B., & Hoffmann, M. 2024, Ann. Appl. Statist., 18, 749

Houk, N., & Smith-Moore, M. 1988, Michigan Catalogue of Two-dimensional Spectral Types for the HD Stars, 4, Declinations -26°.0 to -12°.0.

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Izidoro, A., Raymond, S. N., Morbidelli, A., Hersant, F., & Pierens, A. 2015, ApJ, 800, L22

Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5

Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7

Kipping, D. M. 2013, Mon. Not. R. Astron. Soc. Lett., 434, L51

Kumaraswamy, P. 1980, J. Hydrol., 46, 79

Kurucz, R. 1993, Robert Kurucz CD-ROM, 13

Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83

Leleu, A., Delisle, J. B., Udry, S., et al. 2023, A&A, 669, A117

Li, Y., Brandt, T. D., Brandt, G. M., et al. 2021, AJ, 162, 266

Lo Curto, G., Mayor, M., Benz, W., et al. 2010, A&A, 512, A48

Lo Curto, G., Pepe, F., Avila, G., et al. 2015, The Messenger, 162, 9

Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264

Marmier, M., Ségransan, D., Udry, S., et al. 2013, A&A, 551, A90

Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1

Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619

Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3

Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021, ApJS, 255, 8

Rosenthal, L. J., Knutson, H. A., Chachan, Y., et al. 2022, ApJS, 262, 1

Rutten, R. G. M. 1984, A&A, 130, 353

Santos, N. C., Israelian, G., & Mayor, M. 2001, A&A, 373, 1019

Schlecker, M., Mordasini, C., Emsenhuber, A., et al. 2021, A&A, 656, A71

Skilling, J. 2006, Bayesian Anal., 1, 833

Sneden, C. A. 1973, PhD thesis, University of Texas, Austin

Sousa, S. G. 2014, in Determination of Atmospheric Parameters of B, eds. E. Niemczura, B. Smalley, & W. Pych, 297

Sousa, S. G., Santos, N. C., Israelian, G., Mayor, M., & Monteiro, M. J. P. F. G. 2007, A&A, 469, 783

Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373

Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E., & Israelian, G. 2015, A&A, 577, A67

Sousa, S. G., Adibekyan, V., Delgado-Mena, E., et al. 2021, A&A, 656, A53

Standing, M. R., Triaud, A. H. M. J., Faria, J. P., et al. 2022, MNRAS, 511, 3571

Stevenson, A. T., Haswell, C. A., Faria, J. P., et al. 2025, MNRAS, 539, 727

Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2022, ApJ, 927, 31

Teske, J., Díaz, M. R., Luque, R., et al. 2020, AJ, 160, 96

Tinney, C. G., Butler, R. P., Marcy, G. W., et al. 2002, ApJ, 571, 528

Udry, S., Mayor, M., Naef, D., et al. 2000, A&A, 356, 590

Van Zandt, J., & Petigura, E. A. 2024, AJ, 168, 268

Van Zandt, J., Petigura, E. A., Lubin, J., et al. 2025, arXiv e-prints [arXiv:2501.06342]

Vines, J. I., & Jenkins, J. S. 2022, MNRAS, 513, 2719

Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9

Zhu, W. 2024, Res. A&A, 24, 045013

Zhu, W., & Wu, Y. 2018, AJ, 156, 92
## Appendix A: Extended posteriors table

<table><caption>Table A.1. Maximum likelihood solution and kima posteriors for HD 23079, HD 196067, HD 86226.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">unit</th><th colspan="2">HD 23079</th><th colspan="2">HD 196067</th><th colspan="2">HD 86226</th></tr><tr><th>max-likelihood<sup>a</sup></th><th>kima posterior<sup>b</sup></th><th>max-likelihood<sup>a</sup></th><th>kima posterior<sup>b</sup></th><th>max-likelihood<sup>a</sup></th><th>kima posterior<sup>b</sup></th></tr></thead><tbody><tr><td>v<sub>sys.</sub></td><td>m s<sup>-1</sup></td><td>573.801 ±0.319</td><td>574.0 ± 0.5</td><td>-10945.854 ±0.922</td><td>-10945.0<sup>+1.5</sup><sub>-1.6</sub></td><td>19776.047 ±0.825</td><td>19776.5<sup>+1.6</sup><sub>-1.5</sub></td></tr><tr><td>v<sub>CORALIE98</sub></td><td>m s<sup>-1</sup></td><td>65.81 ±1.30</td><td>65.9 ± 1.4</td><td>64.74 ±2.26</td><td>65.5<sup>+3.0</sup><sub>-2.8</sub></td><td>-31.82 ±3.01</td><td>-33<sup>+4</sup><sub>-5</sub></td></tr><tr><td>v<sub>CORALIE07</sub></td><td>m s<sup>-1</sup></td><td>70.44 ±1.48</td><td>70.6 ± 1.5</td><td>61.83 ±1.35</td><td>62.1<sup>+1.6</sup><sub>-1.5</sub></td><td>-36.89 ±1.31</td><td>-37.6<sup>+1.7</sup><sub>-2.9</sub></td></tr><tr><td>v<sub>CORALIE14</sub></td><td>m s<sup>-1</sup></td><td>95.196 ±0.634</td><td>95.2<sup>+0.7</sup><sub>-0.6</sub></td><td>88.36 ±1.25</td><td>88.7<sup>+1.5</sup><sub>-1.4</sub></td><td>-17.26 ±1.05</td><td>-17.6<sup>+1.3</sup><sub>-1.5</sub></td></tr><tr><td>v<sub>HARPS03</sub></td><td>m s<sup>-1</sup></td><td>-12.112 ±0.525</td><td>-11.9 ± 0.7</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>v<sub>PFS1</sub><sup>c</sup></td><td>m s<sup>-1</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>-160.29 ±1.08</td><td>-160.7<sup>+1.7</sup><sub>-2.2</sub></td></tr><tr><td>v<sub>PFS2</sub><sup>c</sup></td><td>m s<sup>-1</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>-162.46 ±0.832</td><td>-163.3<sup>+1.6</sup><sub>-2.1</sub></td></tr><tr><td>σ<sub>CORALIE98</sub></td><td>m s<sup>-1</sup></td><td>6.46 ±1.10</td><td>6.8<sup>+1.4</sup><sub>-1.2</sub></td><td>6.44 ±1.85</td><td>9.2<sup>+2.6</sup><sub>-2.9</sub></td><td>7.32 ±2.20</td><td>8.3<sup>+3.1</sup><sub>-2.4</sub></td></tr><tr><td>σ<sub>CORALIE07</sub></td><td>m s<sup>-1</sup></td><td>5.28 ±1.21</td><td>5.4<sup>+1.6</sup><sub>-1.3</sub></td><td>8.91 ±1.01</td><td>8.7<sup>+1.1</sup><sub>-1.0</sub></td><td>7.067 ±0.914</td><td>7.1<sup>+1.1</sup><sub>-0.9</sub></td></tr><tr><td>σ<sub>CORALIE14</sub></td><td>m s<sup>-1</sup></td><td>3.019 ±0.651</td><td>3.0<sup>+0.7</sup><sub>-0.8</sub></td><td>4.86 ±1.02</td><td>4.8<sup>+1.3</sup><sub>-1.1</sub></td><td>5.287 ±0.826</td><td>5.1<sup>+1.0</sup><sub>-0.9</sub></td></tr><tr><td>σ<sub>HARPS03</sub></td><td>m s<sup>-1</sup></td><td>1.684 ±0.266</td><td>1.9<sup>+0.4</sup><sub>-0.3</sub></td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>σ<sub>HARPS15</sub></td><td>m s<sup>-1</sup></td><td>1.253 ±0.174</td><td>1.46<sup>+0.27</sup><sub>-0.26</sub></td><td>2.381 ±0.215</td><td>2.48<sup>+0.33</sup><sub>-0.25</sub></td><td>1.435 ±0.163</td><td>1.54<sup>+0.21</sup><sub>-0.18</sub></td></tr><tr><td>σ<sub>PFS1</sub></td><td>m s<sup>-1</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>4.044 ±0.522</td><td>4.1<sup>+0.7</sup><sub>-0.5</sub></td></tr><tr><td>σ<sub>PFS2</sub></td><td>m s<sup>-1</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>1.590 ±0.408</td><td>1.8<sup>+0.6</sup><sub>-0.4</sub></td></tr><tr><td>P<sub>b</sub></td><td>d</td><td>735.740 ±0.220</td><td>735.7<sup>+0.4</sup><sub>-0.4</sub></td><td>3413.7 ±10.5</td><td>3397<sup>+17</sup><sub>-14</sub></td><td>1609.5 ±19.4</td><td>1616<sup>+46</sup><sub>-28</sub></td></tr><tr><td>U<sub>0, b = M<sub>0, b</sub> + ω<sub>b</sub><sup>d</sup></sub></td><td>deg</td><td>107.069 ±0.563</td><td>107.0 ± 0.9</td><td>-50.31 ±1.45</td><td>-53.0<sup>+2.6</sup><sub>-2.2</sub></td><td>-70.55 ±8.07</td><td>-66<sup>+16</sup><sub>-12</sub></td></tr><tr><td>K<sub>b</sub></td><td>m s<sup>-1</sup></td><td>54.386 ±0.248</td><td>54.4<sup>+0.3</sup><sub>-0.4</sub></td><td>82.57 ±1.91</td><td>80<sup>+3</sup><sub>-4</sub></td><td>7.653 ±0.539</td><td>7.8<sup>+0.8</sup><sub>-0.6</sub></td></tr><tr><td>e<sub>b</sub></td><td></td><td>0.07348 ±0.00651</td><td>0.072<sup>+0.012</sup><sub>-0.011</sub></td><td>0.56292 ±0.00703</td><td>0.561<sup>+0.011</sup><sub>-0.012</sub></td><td>0.0389 ±0.0677</td><td>0.07<sup>+0.09</sup><sub>-0.05</sub></td></tr><tr><td>ω<sub>b</sub><sup>e</sup></td><td>deg</td><td>33.20 ±4.92</td><td>29<sup>+15</sup><sub>-7</sub></td><td>142.835 ±0.966</td><td>142.0<sup>+1.5</sup><sub>-1.8</sub></td><td>201.4 ±70.4</td><td>37<sup>+100</sup><sub>-180</sub></td></tr><tr><td>m<sub>b</sub> sin i<sub>b</sub><sup>f</sup></td><td>M<sub>J</sub></td><td>2.66 ±0.05</td><td>2.66 ± 0.05</td><td>5.58 ±0.44</td><td>5.42<sup>+0.20</sup><sub>-0.26</sub></td><td>0.46 ±0.04</td><td>0.49<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>P<sub>c</sub></td><td>d</td><td>5.748913 ±0.000204</td><td>5.74894<sup>+0.00022</sup><sub>-0.00021</sub></td><td>4.601088 ±0.000192</td><td>4.60111 ± 0.00024</td><td>3.984436 ±0.000217</td><td>3.98450<sup>+0.00023</sup><sub>-0.00022</sub></td></tr><tr><td>U<sub>0, c = M<sub>0, c</sub> + ω<sub>c</sub><sup>d</sup></sub></td><td>deg</td><td>155.61 ±9.54</td><td>156 ± 10</td><td>-3.3 ±14.4</td><td>-1<sup>+18</sup><sub>-17</sub></td><td>-17.0 ±22.3</td><td>-10<sup>+24</sup><sub>-22</sub></td></tr><tr><td>K<sub>c</sub></td><td>m s<sup>-1</sup></td><td>2.702 ±0.252</td><td>2.65<sup>+0.26</sup><sub>-0.27</sub></td><td>3.762 ±0.396</td><td>3.6 ± 0.4</td><td>2.866 ±0.253</td><td>2.82<sup>+0.25</sup><sub>-0.26</sub></td></tr><tr><td>e<sub>c</sub></td><td></td><td>0.097 ±0.124</td><td>0.07<sup>+0.09</sup><sub>-0.05</sub></td><td>0.235 ±0.118</td><td>0.12<sup>+0.12</sup><sub>-0.09</sub></td><td>0.0619 ±0.0817</td><td>0.05<sup>+0.07</sup><sub>-0.04</sub></td></tr><tr><td>ω<sub>c</sub><sup>e</sup></td><td>deg</td><td>283.4 ±47.2</td><td>247<sup>+64</sup><sub>-166</sub></td><td>295.5 ±23.8</td><td>278<sup>+41</sup><sub>-104</sub></td><td>144.6 ±91.4</td><td>170<sup>+108</sup><sub>-94</sub></td></tr><tr><td>m<sub>c</sub> sin i<sub>c</sub><sup>f</sup></td><td>M<sub>⊕</sub></td><td>8.26 ±0.79</td><td>8.1 ± 0.8</td><td>10.4 ±1.39</td><td>10.1 ± 1.2</td><td>7.45 ±0.7</td><td>7.7 ± 0.7</td></tr><tr><td>reference epoch</td><td></td><td>2 455 000 BJD</td><td></td><td>2 455 000 BJD</td><td></td><td>2 455 000 BJD</td><td></td></tr><tr><td>N</td><td></td><td>199</td><td></td><td>187</td><td></td><td>246</td><td></td></tr><tr><td>log L<sub>max</sub></td><td></td><td>-537.8944</td><td></td><td>-569.7350</td><td></td><td>-690.3457</td><td></td></tr><tr><td>log Z</td><td></td><td>-625.97</td><td></td><td>-644.7</td><td></td><td>-767.46</td><td></td></tr><tr><td>I</td><td></td><td>77.19</td><td></td><td>67.56</td><td></td><td>70.37</td><td></td></tr></tbody></table>

Notes. (a) The provided values are the maximum-likelihood estimates, and the errorbars are the 1-σ uncertainties estimated from the second derivatives of the log-likelihood. (b) The provided values are the posterior medians and the errorbars correspond to the 16% and 84% quantiles.

(c) An arbitrary offset of 19619.947 m/s (the systemic RV of HD 86226 measured by Gaia) was added to the PFS RVs, since the data published by Teske et al. (2020) were already mean-subtracted. (d) Mean argument of latitude at the reference epoch (U₀ = M₀ + ω = λ₀ - Ω, where M₀, λ₀, ω, and Ω are the mean anomaly, mean longitude, argument of periastron, and longitude of the ascending node, respectively). (e) Argument of periastron of the stellar orbit. (f) The uncertainties on the planetary masses take into account the stellar mass formal uncertainties, but not stellar models systematic errors (see Table 1).
**Appendix B: Compatibility limits without PFS data**

The compatibility limits shown in Fig. 6 (right) for HD 86226 include PFS radial velocities in addition to CORALIE and HARPS data. This is useful to assess whether additional planets could still hide in this system. However, in order to assess the sensitivity of our observing programme to detect small planets, we additionally produced a compatibility map using only CORALIE and HARPS data, which is the typical case for most targets in our programme.

<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Compatibility limits of HD 86226 including (top, see also Fig. 6) or not (bottom) PFS RVs in the analysis.</figcaption></figure>

Both compatibility maps are shown in Fig. B.1 for comparison. Adding the PFS data improves the long-period planets detectability but only marginally the short periods planets detectability. Short periods are well covered by HARPS measurements which have a significantly better precision than PFS measurements (see Table 2). Thus the addition of PFS data does not improve much the SNR at short period. For longer periods, including PFS measurements allows to have a longer baseline covered with high-precision RV instruments.

**Appendix C: Inner planets transit-search with TESS**

In this appendix we present the search for potential transits of the inner planets detected for each star. For this analysis, we fix the stellar parameters at their median value reported in Table 1. For each star, the raw PDCSAP flux was downloaded using the `lightkurve`<fnref n="5" /> package. We obtained sectors 3, 4, 30 and 31 for

HD 23079, sectors 13, 27, 66 and 67 for HD 196067, and sectors 9, 35 and 62 for HD 86226. Then, we split the time series into half-sectors in order to model the noise on each independently. We detrend the data with a B-spline using the `linmarg`<fnref n="6" /> Python package (see Leleu et al. 2023). The spacing of the nodes of the B-spline is chosen to be 3 times the transit duration of a planet on a circular orbit at the period of the RV detection and with a null impact parameter.

We first performed a BLS search around the period found in RV with a $3\sigma_{P_c}$ window in period ($\sigma$ is chosen as the max of the two errors reported in Table A.1) around the median $\bar{P}_c$, to see if there is an obvious transit in the data. We found HD 86226 c with properties similar to the one reported in Teske et al. (2020), but no signal with SNR more than 4 in this pre-detrended dataset for HD 23079 c and HD 196067 c.

<figure><img src="image_11.png" /><figcaption>**Fig. C.1.** Model comparison map for HD 23079 c, for radius ranging from 0.5 to 3.5 $R_⊕$ and impact parameters ranging from 0.1 to 0.9. The colors indicate the value of the $\Delta\text{BIC}$, red indicating that the null model is favoured, while blue indicates that the planet model is favoured. a value of $|\Delta\text{BIC}| > 20$ is usually considered to be a compelling evidence for the model.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. C.2.** Model comparison map for HD 196067 c, see Fig. C.1 for the text for more detail.</figcaption></figure>

<footnote n="5">⁵ https://docs.lightkurve.org/</footnote>

<footnote n="6">⁶ https://gitlab.unige.ch/delisle/linmarg</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. C.3.** Model comparison map for HD 86226 c, see Fig. C.1 for a text

for more detail.</figcaption></figure>

We then computed a model comparison map for each planet
with a simultaneous noise and planet modelling. To do so, we
first identify potentially problematic parts of the light curves
using the pre-detrended data. We start by flagging the data points
that are at a larger separation from the mean value than 4 times
the standard deviation. We then binned the data with a 30 min
cadence. We then checked, for each half-sector, if the standard
deviation of the flux was 2.5 larger than the minimum standard
deviation of all the binned half-sector for that star. This notably
lead to the removal of sector 66 for HD 196067 which had a
much larger scatter than the others. All other half-sectors were
kept. We then identified all time stamps that were removed in
the detrended light curves, and removed it from the raw data.
Restarting from the masked raw data, we then estimated the like-
lihood of the null model L₀, only applying the detrending. Then,
for a grid of planetary radii R_c and impact parameters b_c, we
explored a grid of orbits ranging the 3σ interval both in period
and λ₀ for each inner planet. The step in λ₀ is set such that the
time of transit around the reference epoch varies by 10 min steps.
The step in period is set such as the last transit varies by steps
of 15 min (δP = 6σp_c (t_max - t_min)/p_c)/96. Then, for each cou-
ple (R_c, b_c), we save the highest likelihood L_p in the grid (t_0, c,
P_c), with a simultaneous modelling of the transits of the planet
and the noise model (the aforementioned B-splines). We then
compute the ΔBIC = 2(L_p - L₀) - 4 log(N_bins).

Figures C.1, C.2 and C.3 show the model comparison maps for HD 23079 c, HD 196067 c and HD 86226 c, respectively. For HD 23079 c and HD 86226 c, the interpretation of the maps are straightforward: existing data exclude a planetary transit for HD 23079 c, while the known ~2R<sub>⊕</sub> planet is favoured for HD 86226 c (Teske et al. 2020). For HD 196067 c, we notice a ΔBIC ~ 48 for R<sub>c</sub> ~ 1.9R<sub>⊕</sub> and an impact parameter of 0.9. We investigated this solution by estimating ΔBIC per half-sector instead of the whole dataset. It turns out that ΔBIC is only positive for the first half of sector 13, and negative for the 5 other half-sectors. Upon visual inspection of the lightcurve, a single dip of the flux over a timescale of ~1h seems to be responsible for the positive value of the ΔBIC in the first half of sector 13. We therefore conclude that the planet is not transiting.