© The Authors 2026

1 Introduction

Metal-rich stars have played a central role in exoplanet discovery since the earliest detections of gas giants, and the planet-metallicity correlation remains one of the most robust demographic trends: gas giants are markedly more common around metal-rich FGK stars (Gonzalez 1997; Santos et al. 2004; Fischer & Valenti 2005; Winn 2018). Today, there are over 6000 confirmed exoplanets, the vast majority found by radial velocities (RVs) and transit photometry. However, the sensitivity of both techniques declines with orbital separation a (Zhu & Dong 2021): the geometric transit probability scales as R_*/a, and long-period RV signals require many years of high-precision monitoring to resolve slowly varying signals with a low amplitude (Cumming et al. 2008). As a result, the population of cold-orbiting giants beyond a few AU remains comparatively unexplored.

When we refer to exoplanets here, we adopt the term ‘cold’ for those with a semi-major axis of 1-10 AU, ‘wide-orbiting’ for >10 AU, and ‘Jupiter analogue’ for 3-7 AU and a mass of 0.3-13 M_J (Bryan et al. 2019). Recent occurrence-rate analyses showed that Jupiters with orbital periods >100 d occur around Sun-like stars at 6-7%. They are more common by nearly an order of magnitude than hot Jupiters (Wittenmyer et al. 2020). Moreover, the relative occurrence of hot to cold Jupiters appears to decrease with stellar mass, consistent with expectations from classical core accretion and migration timescales (Gan et al. 2024). These trends highlight the need for long-baseline surveys capable of detecting and characterising cold giants.

While precise RVs constrain orbital parameters and minimum masses M_pl sin i, they do not measure inclination, and thus, they cannot yield true masses by themselves. Absolute astrometry, on the other hand, is directly sensitive to the sky-projected stellar reflex motion and provides dynamical masses and inclinations, albeit with sensitivity that degrades with distance d. In particular, the signal-to-noise ratio (S/N) scales as $\mathrm{S}/{\mathrm{N}}_{\mathrm{A}\mathrm{M}}\propto \frac{M_{\mathrm{p}\mathrm{l}}\cdot a}{M_{\star }\cdot d}$, while for RVs, this ratio scales as $\mathrm{S}/{\mathrm{N}}_{\mathrm{R}\mathrm{V}}\propto \frac{M_{\mathrm{p}\mathrm{l}}}{\sqrt{a\cdot M_{\star }}}$.

Combining RVs with absolute astrometry breaks key degeneracies, enhances the S/N, and enables a robust characterisation of long-period companions even when a single technique samples only a fraction of the orbit. With HIPPARCOS (Hipp) and Gaia providing high-precision positions, proper motions, and accelerations over a baseline of −24 yr, the RV+astrometry synergy has matured to a powerful pathway for measuring true masses of cold giants and identifying a precise dynamical mass for Jupiter analogues (Perryman et al. 2014; Kervella et al. 2022; Feng et al. 2019; Brandt 2021).

A paradigm shift is imminent within the precision-RV community: long-period RV candidates will increasingly be treated as testable astrometric hypotheses and systematically checked for consistency with Gaia. In practice, an RV posterior predicts an on-sky reflex signal (amplitude and orientation) that Gaia can confirm, refute, or tighten into a true-mass measurement. Recent forecasts suggest that Gaia astrometry may yield 7500 ± 2100 planet detections in DR4 and 120 000 ± 22 000 in DR5, with −1900 ± 540 (DR4) and 38 000 ± 7300 (DR5) having masses and periods measured to be better than 20%. Most of these detections are expected to be super-Jupiters on −2-5 AU orbits around nearby GKM stars (Lammers & Winn 2026), overlapping with the parameter space targeted by long-baseline Doppler surveys.

The connection between outer cold giants and small inner planets is crucial for understanding system architectures. Early observational work reported that systems with inner super-Earths frequently host cold giants (roughly one in three), implying a positive correlation (Zhu et al. 2018; Bryan et al. 2019). Subsequent studies found lower conditional rates or no excess in carefully vetted samples (Bonomo et al. 2023). A more nuanced, mass-dependent picture is emerging: new analyses indicate that Saturn-mass cold giants are positively correlated with inner superEarths, whilst the most massive super-Jupiters are not, potentially reconciling previous discrepancies and informing models of accretion efficiency and dynamical histories in multi-planet systems (Lefèvre-Forjàn & Mulders 2025). Establishing how common true Jupiter analogues are and measuring their dynamical masses directly bears on how typical the Solar System is among planetary systems.

Cold giants regulate the assembly and survival of inner planetary systems, sculpt debris belts, and provide boundary conditions for planet formation theories (migration histories, core accretion Ida et al. 2013; Mordasini et al. 2012). Measuring their occurrence, mass function, and orbital architectures tests these theories where they differ most, informs the demographics accessible to direct imaging, and provides prime targets for follow-up. The newly developed RV-astrometry coupling is explicitly built into future facilities: the Tianyu transit project plans confirmation and characterisation of long-period giants via Gaia astrometry and RV follow-up, leveraging multi-survey baselines (Feng et al. 2024b); and Habitable Worlds Observatory (HWO; Dressing et al. 2024) studies emphasise that precursor RVs together with a modest number of astrometric epochs can deliver full orbital solutions for cold giants and habitable-zone planets, enabling efficient direct-imaging campaigns (Sagynbayeva et al. 2025).

The Chile-Hertfordshire ExoPlanet Survey (CHEPS) was initiated in 2009 to target inactive, metal-rich stars with high-precision RVs from HARPS and CORALIE, explicitly leveraging the planet-metallicity correlation to find gas giants efficiently (Jenkins et al. 2008, 2009, 2017). By extending the CHEPS observations through 2025, we assembled a 16-year time baseline (2009-2025) ideally suited to unveil long-period curvatures and complete (or near complete) orbits for cold exoplanets. In this paper, we leverage this long CHEPS baseline with modern astrometric constraints. We (1) identify five CHEPS targets exhibiting long-term RV curvature (HIP 8923, HIP10090, HIP21850, HIP39330, and HIP98599); (2) model each system with a joint RV+astrometry framework using EMPEROR (Peña R & Jenkins 2026); and (3) determine the true masses and orbital parameters of the suspected cold giants. Beyond target-level results, our analysis serves as an end-to-end validation for the EMPEROR approach for combined RV+astrometry inference on long period planets. The paper is organised as follows: Sect. 2 outlines the joint modelling framework; Sect. 3 summarises datasets and stellar parameters; Sect. 4 presents system-by-system fits; Sect. 5 discusses and quantifies the benefits of the joint modelling, places the discovered systems in the context of the overall exoplane-tary census, and addresses possible drawbacks in this study; and Sect. 6 summarises the key points of this research.

2 Overview of the modelling framework

EMPEROR (Exoplanet Mcmc Parallel tEmpering for Rv Orbit Retrieval; Peña R & Jenkins 2026) is a modular Python-based pipeline that searches RV time series for planetary signals and performs Bayesian model comparison to select the preferred solution. Posteriors are sampled with reddemcee (Peña R & Jenkins 2026), an adaptive parallel tempering (APT) MCMC sampler (Earl & Deem 2005; Vousden et al. 2016) designed to explore broad multi-modal posteriors efficiently and mitigate local trapping, significantly accelerating convergence rates while improving the reliability of global inferences. The framework estimates the marginalised likelihood (also known as evidence) through a hybrid approach that combines thermodynamic integration (Gelman & Meng 1998) and stepping stones (SS, Xie et al. 2011) with realistic uncertainties. EMPEROR supports multi-instrument, multi-planet analyses, and flexible noise models including (when needed) stellar activity correlations and Gaussian processes. In what follows, we summarise the Bayesian formulation and models used, and describe how we incorporated the HIPPARCOS-Gaia differencing framework in the joint likelihood (Feng et al. 2019).

2.1 Bayesian formulation

The Bayes theorem defines the posterior p(θ | D, M) of parameters θ, given data D and model M as

$$p(\theta \theta \mid D,M)=\frac{p(\theta \theta \mid M)\cdot p(D\mid \theta \theta ,M)}{p(D\mid M)},$$(1)

where p(θ | M) is the prior, p(D | θ, M) the likelihood, and p(D|M) the Bayesian evidence normalising the posterior (Robert 2007; Fong & Holmes 2019; Lotfi et al. 2022). The evidence enables objective model comparison via Bayes’ factors. In our case, the likelihood is joint, coupling RVs and astrometry (subscript AM) through the shared orbital parameters θ,

$$p(D\mid \theta \theta ,M)=p(D_{\mathrm{R}\mathrm{V}}\mid \theta \theta ,M)\times p(D_{\mathrm{A}\mathrm{M}}\mid \theta \theta ,M).$$(2)

2.2 Radial velocity model

The deterministic RV signal is a sum of Keplerian curves,

$$\mathcal{K}(t)=\sum _{j=1}^{N_{\mathrm{p}\mathrm{l}}}{K}_j\cdot [\cos ({\nu }_j(t)+{\omega }_j)+e_j\cos ({\omega }_j)],$$(3)

where each planet j has a period P_j, semi-amplitude K_j, eccentricity e_j, longitude of periastron ω_j, time of periastron passage T_0,j and true anomaly ν_j. Instrumental zero-points γ_0,INS, long-term trends γ, and (optionally) linear correlations with stellar activity indices C_A,INS ∙ A_i,INS (with A_i,INS denoting the measured activity index and C_A,INS the coefficient) are included additively into Eq. 3. Two orbital angles remain unconstrained by RVs alone, the inclination i_j (tilt with respect to the plane of sky) and the longitude of ascending node Ω_j (plane-of-sky orientation). Consequently, RVs measure only the minimum mass M_pl,j sin i_j via the mass function

$$\frac{P_j{K}_j^3\sqrt{(1-e_j^2{)}^3}}{2\pi G}=\frac{(M_{\mathrm{p}\mathrm{l},j}\cdot \sin i_j{)}^3}{(M_{\mathrm{p}\mathrm{l},j}+M_{\star }{)}^2}.$$(4)

2.3 Astrometric model (Hipparcos-Gaia differencing with IAD+GOST)

We adopted the astrometric differencing approach of Feng et al. (2019, 2023, 2024a). Rather than relying only on catalogue-level proper motion differences, we (i) fitted the HIPPARCOS intermediate astrometric data (IAD, the along-scan measurements at ∼ 1990-1993; van Leeuwen 2007); and (ii) propagated and projected our model at the Gaia reference epoch using scan angles from the Gaia observation forecast tool (GOST; ESA/DPAC 2016a,b). This IAD+GOST approach preserves the HIPPARCOS epoch information and properly matches Gaia ’s along-scan sensitivity, avoiding the information loss that can occur when collapsing to catalogue-level proper motion anomalies.

We propagated a five-parameter astrometric model between HIPPARCOS and Gaia epochs (∆α, ∆δ, ∆μ_α, ∆μ_δ, ∆ϖ), including perspective terms from catalogue RV and parallax. We added the star’s 3D reflex motion due to the companion(s) and projected it onto the sky. Internally, we used Thiele-Innes constants to map Keplerian elements to sky-plane displacements (for a detailed description, see Xiao et al. 2024). For every HIPPARCOS scan, we computed the predicted along-scan displacement from the model (barycentre + reflex) at the scan time and angle and evaluated a Gaussian likelihood using the IAD uncertainties plus an additive jitter term J_hipp. On the other hand, at the Gaia reference epoch, we projected the model into Gaia’s along-scan direction using the GOST angles and parallax factors and formed the position and proper-motion differences relative to the HIPPARCOS solution linearly propagated to the Gaia epoch. The astrometric likelihood corresponds to the product of the HIPPARCOS-IAD term and the Gaia difference term. We included per-catalogue offsets and small frame-rotation jitters.

After linear-space motion was removed, a companion induced a sky-plane acceleration g that left two signatures over the HIPPARCOS-Gaia epoch separation ∆t,

$$\mathrm{\Delta }r\approx \frac{1}{\text{2}}g(\mathrm{\Delta }t{)}^2,\mathrm{\Delta }\mu \approx g\mathrm{\Delta }t,$$(5)

so jointly modelling positional (∆r) and proper-motion (∆μ) offsets is especially sensitive to long-period companions. Using the actual scan geometries (GOST) reduces susceptibility to scanangle-dependent systematics (e.g., Holl et al. 2023). Then, in this framework, the GOST model corresponds to the geometric projection layer that converts the orbital model in sky coordinates into Gaia ’s effective along-scan observable using GOST scan geometry, while the GDR2/3 model corresponds to the statistical likelihood layer that compares that predicted mean vector to the catalogue-level Gaia solution using the reported covariance matrix. Using multi-epoch Gaia astrometry (GDR2 and 3) can help us to differentiate between prograde and retrograde orbits, breaking the inclination degeneracy (Xiao et al. 2024).

2.4 Priors

Our prior parameters were purposely flat for most RV and astrometric parameters. For long-period giants, the empirical eccentricity distribution is broad (Kipping 2013; Stevenson et al. 2025), but to regularise poorly constrained cases and to avoid bias towards spuriously high e in low-SNR regimes, we used a truncated normal prior ~N(0,0.3) ∈ [0,1], penalising extremely higher eccentricities, but allowing them if supported by the data (see Tuomi et al. 2013). RV jitters use a wide truncated normal centred at zero, with the scale set by the average precision of the instrument, at either σ = 5 or 10 m s⁻¹. HIPPARCOS additive astrometric jitter J_hipp~N(1, 3) and fractional Gaia jitter S_gaia~N(1,0.1) absorb modest under- or over-dispersions.

For systems with two planets, we re-parametrised inclinations as (I₁,I₂) → (I₀,I_∆), with

$$I_1=I_0-I_{\mathrm{\Delta }},I_2=I_0+I_{\mathrm{\Delta }}.$$(6)

We favoured quasi-coplanarity via I_Δ~N(0,5.73°), which is consistent with multi-planet mutual inclination dispersion (Fabrycky & Winn 2009; Zhu et al. 2018; He et al. 2020), while allowing separated inclinations if supported by the data. The reference inclination I₀~I(0,180°) uses an isotropic prior, uniform in cos I₀ over [0,1].

3 Target selection and observations

Since its inception, CHEPS has targeted quiet, metal-rich FGK stars (Jenkins et al. 2008; Jenkins & Jordán 2011) with the goal of discovering giant planets (Jenkins et al. 2009, 2017). For the present study, we selected CHEPS stars that showed significant long-term RV curvature or acceleration in the 2009-2021 time series and continued monitoring them through 2025, thereby extending individual baselines up to ∼16 years. In addition to new CHEPS observations, we assembled all available high-precision archival RVs for these stars from HARPS, CORALIE, UCLES, and HIRES. For this first batch of CHEPS systems, we selected inactive stars (log R′_HK ≤ –4.75), with low S_MW values (S_MW ≤ 0.201) that were also slow rotators (v sin i ≤ 5 km s⁻¹). Furthermore, we considered targets with a total RV baseline over 14 years and more than 16 CORALIE measurements (see the full stellar parameters in Table 1).

Our RV sample combines high-precision spectra from the High Accuracy Radial Velocity Planet Searcher (HARPS; Pepe et al. 2000) on the ESO 3.6 m telescope at La Silla Observatory; CORALIE (Udry et al. 2000) on the 1.2 m Euler telescope at La Silla Observatory; University College London Echelle Spectrograph (UCLES; Diego et al. 1990) at the Anglo-Australian Telescope (AAT); and the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) at the Keck Observatory. The HARPS data were retrieved from the public HARPS RVBank (Trifonov 2019; Perdelwitz et al. 2024), whilst UCLES and HIRES archival velocities were gathered from published releases and institutional archives.

HARPS experienced a fibre upgrade (2 June 2015) and a warm-up (23 March 2020) post-COVID shutdown (Lo Curto et al. 2015). Therefore, we treated HARPS with three different offsets: H (pre-2015), H15 (2015-2020), and H20 (post-2020).

CORALIE underwent a major upgrade in 2014. The circular fibre link was replaced with octagonal fibres, and a Fabry-Pérot calibration unit was added (Rickman et al. 2019). We adopted three subsets consistent with CHEPS processing: COR07 (pre-upgrade), COR14T (post-upgrades, reduced with TERRA; Anglada-Escudé & Butler 2012), and COR14 (post-upgrades, reduced with the CORALIE DRS).

We used the HIPPARCOS IAD from the new 2007 reduction (Perryman et al. 1997; van Leeuwen 2007), obtaining the data through htof (Brandt 2021), and Gaia DR3 five-parameter solutions (Gaia Collaboration 2023). To project our orbital models into Gaia’s measurement geometry at the Gaia reference epoch, we queried GOST for scan angles. A summary table with the number of measurements per target can be found in Table C.1.

4 Analysis and results

Within the EMPEROR framework, joint astrometry is straightforward to use, and a step-by-step tutorial can be found in the online documentation¹.

For each target, we examined the RVs with generalised Lomb-Scargle periodograms (GLS; Lomb 1976; Scargle 1982), as well as with correlograms, calculating the Pearson-correlation coefficients between the RVs and stellar activity indices. We then began to fit models with EMPEROR, first by testing the null hypothesis that comprises just a mean and white noise (hereinafter 0K). After this, we incrementally added Keplerian signals (models dubbed 1K, 2K, and so on) and compared the Bayesian evidence ln Z with the previous model; with a difference of ∆ln Z ≥ 5, we selected the new model. Each model was run seven times to ensure consistent results. The values we present for the model comparison and parameter estimation correspond to the run with the median evidence. The hyper-parameter selection for our runs was set to 20 temperatures, 256 walkers, 80 000 sweeps, and 1 step for a grand total of 409 600 000 samples. Ladder adaptation was made with a uniform swap acceptance rate kernel, with a rate ν₀ = 1.0 and a timescale τ₀ = 100. We left the adaptation for 40 000 sweeps and then froze the ladder. This initial phase we considered as burn-in, so these samples were dropped before any statistical inferences were made, preserving detailed balance so that standard ergodicity theorems apply.

As a validation case, we first analysed HIP 21850, a well-studied system hosting two confirmed exoplanets. We took as references the parameter solutions of Wittenmyer et al. (2017b, RV-only) and Feng et al. (2022, joint RV+astrometry) (W17 and F22, respectively). A summary of planet signatures for all systems can be found in Table 2.

4.1 HIP 21850

The RVs of planet HD 30177 present a highly significant peak in the GLS periodogram at above the 99.9% false-alarm probability (FAP) line at 2629 d. The GLS of the full width at half-maximum (FWHM) has a marginally significant peak between the 90% and 99.0% FAP lines at 2596 d (see Fig. 1), while the bisector inverse slope (BIS) presents a highly significant peak at 2671 d. The window function presents a dominant peak at 2960 d. On the other hand, the correlogram reveals weak to moderate correlations (|ρ| ∈ [0.20,0.60]) between RVs and indices: ρ_RV,FW=-0.31 for FWHM and ρ_RV,BIS=0.45 for BIS (see Fig. 2).

With the EMPEROR framework, we started by analysing RVs alone. The single-Keplerian solution (1K) finds an orbital period of $P_1={2534}_{-11}^{+22}$ d, and adding a second Keplerian (2K) yields periods $P_1={2539.8}_{-4.7}^{+4.1}$ d and $P_2=15{900}_{-4800}^{+2500}$ d. The model comparison is decisive for the 2K model, as the Bayes factor (BF) $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=102.4$ strongly favours it, with a Bayesian evidence of $\ln \hat{\mathcal{Z}}{\mid }_{2\mathcal{K}}=-389.5\pm 0.3$ over the 1K model's $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-491.9\pm 0.2$. Each planet has minimum mass and semi-major axis estimates of m sin i₁ = ${8.2}_{-0.5}^{+0.1}{M}_{\mathrm{J}}$, m sin i₂ = ${6.6}_{-2.5}^{+1.4}$ M_J, a₁ = ${3.65}_{-0.11}^{+0.02}$, and a₂ = ${12.5}_{-2.7}^{+1.4}$.

The preference is the same for EMPEROR with RV+astrometry, $\ln \hat{\mathcal{Z}}{\mid }_{2\mathcal{K}}=-713.4\pm 0.3$ over $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-817.5\pm 0.3$, with a BF of $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=104.10$ (we extensively discuss the BF differences with the addition of astrometry in Sect. 5.1). This time, we find $P_1={2539}_{-3}^{+1}$ d and $P_2=11{320}_{-940}^{+490}$ d, with true masses and semi-major axes $M_1={8.25}_{-0.85}^{+1.05}{M}_{\mathrm{J}}$, $M_2={4.67}_{-0.43}^{+0.44}{M}_{\mathrm{J}}$, $a_1={3.49}_{-0.14}^{+0.18}$ and $a_2={9.90}_{-0.62}^{+0.29}$ (see Table 3).

The inclusion of astrometry significantly helps us to constrain the second period, $P_2=15{900}_{-4800}^{+2500}\to 11{320}_{-940}^{+490}$ d, while determining the true masses. The posteriors for the 2K RV-only model reveal that P₂ sits uncomfortably in a posterior with a flattened top between ∼9000-17 000 d, as attested by the loose uncertainties. The RV+astrometry posteriors round significantly, still with a flattened top around 11 000-12 000 d, which adds a clear upper boundary on the solution. Our solution remains fairly similar to those of W17 and F22, considering that with the addition of the CHEPS RV data, we extended the baseline from ∼6000 d to 9070 d. The RV part of the joint-model fit is shown in Fig. 3. The astrometry is shown in Fig. 4. We further discuss the parameter estimates in Sect. 5.2.

For consistency, we again ran the RV residuals through the GLS: the RV field was flattened, with a single marginally significant peak at 90% in the FAP line at 91.1 d. Similarly, the residual correlogram was flat, with ρ_RV,FW=0.01 and ρ_RV,BIS =-0.07. The flat residual periodogram and correlogram, the strong Bayesian preference, and the internal consistency of the P₁ solution are good indicators for a planetary companion: if stellar activity were the origin of the main RV peaks, we would expect coherent residual power and/or correlations. Furthermore, all RV signals we discussed have larger amplitudes than are expected to be produced by stellar activity for these stars (up to a few ms⁻¹; Makarov et al. 2009; Meunier et al. 2010). We further discuss the stellar activity and modelling with GPs with the EMPEROR framework in Peña R & Jenkins (2026).

Finally, is worth discussing that the window function, initially a plausible source for some GLS structure given the peaks at 2960 and 2360 d, is unlikely to shape this solution: for 1K, the fitted P₁ = 2539 d signal sits between the two window function spikes and survives in between when moving onto 2K, and even when astrometry is included in the fit.

Fig. 1 HIP 21850 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and the window function. The FAP lines are included for 0.1, 1, and 10% and are plotted as dashed red, dotted purple, and dotted blue lines, respectively. The circles show the periods with the greatest power, coloured by FAP region. The coloured regions correspond to $P_1={2539}_{-3.1}^{+0.8}$ d (green) and $P_2={11320}_{-940}^{+490}$ (orange).

Fig. 2 HIP 21850 correlogram. Each panel displays the pairwise relation between RVs, FWHM, CCF BIS, and RVe. The Pearson correlation coefficient is denoted by ρ, and the linear trend is shown as a black line.

Fig. 3 HIP 21850 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.

Fig. 4 HIP 21850 best-fit astrometry for the 2K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

4.2 HIP 8923

The planet HIP 8923 is also known as HD 11731. The GLS periodogram (see Fig. A.1) shows three highly significant peaks (>99.9% FAP line) in the RVs, at 5020 d and 314 d, which correspond to the yearly alias between 5020 d and 337 d in the window function, and at 29.7 d (the moon cycle), which is also present in the window function. The correlogram, on the other hand, shows a strong correlation (here defined as |ρ| ∈ [0.60,0.80]) in ρ_RV,FW=-0.63. This correlation might be an artefact due to the uncertainty difference between datasets, since there is a very strong correlation (ρ ∈ [0.80,1]) between the RV uncertainty and the FWHM ρ_RVe,FW=0.97.

Our RV-only search yielded a single planet with an orbital period $P_1={5120}_{-260}^{+460}$ d, with evidence $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-87.64\pm 0.11$, and a BF of $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=16.77$ over the white-noise-only model $\ln \hat{\mathcal{Z}}{\mid }_{0\mathcal{K}}=-104.41\pm 0.09$. By adding the astrometry data, we found $P_1={5160}_{-240}^{+150}$ d, which is consistent with the RV-only model. In this scenario, the evidence is $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-297.97\pm 0.19$, with a $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=75.35$ over the 0K model $\ln \hat{\mathcal{Z}}{\mid }_{0\mathcal{K}}=-373.32\pm 0.23$.

This case strongly demonstrates the value of astrometry when only a part of a long-period orbit is sampled by RVs (see Fig. 5). The relative BF between 0K and 1K increases significantly when astrometry is added, $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=16.77\to 75.35$. Furthermore, the inclusion of the astrometric data shrank the period uncertainty almost by half $P_1={5120}_{-260}^{+460}\mapsto {5160}_{-240}^{+150}$ d, where the relative uncertainty corresponds to 13.97%, and 7.57%, respectively. Similarly, the minimum mass m sin i=${4.30}_{-0.75}^{+0.07}$ M_J becomes an absolute mass of $M_1={9.98}_{-0.16}^{+0.78}{M}_{\mathrm{J}}$, and semi-major axis $a_1={5.89}_{-0.26}^{+0.30}\mapsto {5.90}_{-0.22}^{+0.10}$. The complete astrometric fit can be seen Fig. 6.

Finally, we investigated the peak that appears in the window function at 4818 d. Even though it is ranked 12th in height, we analysed the residuals of the model: the GLS leaves no significant peaks above the 90% FAP line, which is a strong argument that the 5000-day power is well explained by the orbit rather than harmonics or aliases. On the other hand, the residual correlogram shows ρ_res,FW=0.28, a significant shift from the initial ρ_RV,FW=-0.63. Even though this correlation is most likely an artefact due RV uncertainties between datasets, the strong anticorrelation disappears when the Keplerian is removed, which is consistent with the planet driving most of the long-term RV variation, and does not indicate stellar activity. We also note that the yearly alias producing the second highest peak in the RVs at 314 d completely disappears in the residuals.

Fig. 5 HIP 8923 RVs phase-folded at $P={5160}_{-240}^{+150}$ d as circles coloured per instrument with the 1K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right of each plot, we show RV distribution histograms.

Fig. 6 HIP 8923 best fit for the 1K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

4.3 HIP 10090

For HIP 10090, also known as HD 13350, our GLS periodogram reveals a dominant highly significant peak at 2892 d, another at 324 d, and a family of peaks around 25-32 d, corresponding to the lunar cycle. The stellar activity indices show no peaks above the 90% FAP line, and the window function shows peaks at 364 and 29.7 d, corresponding to the yearly and monthly aliases. The correlogram shows no significative correlations with stellar activity (|ρ| < 0.2).

With EMPEROR, the RV-only solution settles at 2K with periods $P_1={2980}_{-100}^{+140}$ d and $P_2={321.3}_{-0.4}^{+0.5}$ d, with evidence $\ln \hat{\mathcal{Z}}{\mid }_{2\mathcal{K}}=-203.84\pm 0.18$, and a BF of 10.92 over the $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}$. The runs including astrometry data settle at 2K as well, with $P_1={2960}_{-100}^{+120}$ d and $P_2={321.8}_{-0.6}^{+0.3}$ d (see Fig. 7), consistent with the RV-only solution. This run has $\ln \hat{\mathcal{Z}}{\mid }_{2\mathcal{K}}=-548.57\pm 0.19$, with a BF of 27.22 over $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}$. By analysing the residuals of this model, the GLS stays flat, without peaks over the 90% FAP line, while the correlations remain insignificant. Once again, astrometry drives to a precise mass estimation m sin i₁ = ${1.53}_{-0.14}^{+0.12}\mapsto M_1={3.87}_{-0.60}{M}_{\mathrm{J}}$, and m sin i₂ = ${0.39}_{-0.03}^{+0.02}\mapsto M_2={0.85}_{-0.12}^{+0.03}{M}_{\mathrm{J}}$.

The relatively high proper motion of HIPPARCOS (red ellipse in Fig. 8) makes this part of the model relatively uninformative compared to Gaia GDR2 and 3. We further discuss this point in Sect. 5.1.

Fig. 7 HIP 10090 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right of each plot, we show RV distribution histograms.

Fig. 8 HIP 10090 best fit for the 2K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

4.4 HIP 39330

For HIP 39330, also known as HD 66653, the GLS periodogram shows a single dominant peak in the RVs (see Fig. A.4) at 4547 d. The FWHM periodogram has significant peaks at 17.9 d and 29.0 d, the latter corresponding to the moon cycle. The correlogram does not reveal significant correlations (|ρ| < 0.2) between RVs and stellar activity (see Fig. B.4).

The RV-only solution settles at 1K with $P_1={4710}_{-150}^{+20}$ d, with $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-202.41\pm 0.14$. The 2K solution presents higher evidence $\ln \hat{\mathcal{Z}}{\mid }_{2\mathcal{K}}=-198.48\pm 0.17$, but it does not comply with our previous requirement of $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=3.93<5$, finding orbital periods $P_1={4640}_{-200}^{+300}$ d, and $P_2={257}_{-75}^{+160}$ d.

The joint RV+astrometry analysis settles at 2K with $\ln \hat{\mathcal{Z}}{\mid }_2=-548.83\pm 0.23$, with a $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=6.97$ over the 1K, with $P_1={4596.9}_{-154.0}^{+1.7}$ d, and $P_2={257.1}_{-1.4}^{+1.6}$ d (see Fig. 9). This system alone of our sample presents equally likely multiple peaks (see Fig. C.1). Our sampler, or any MCMC sampler, is characterised by not returning point estimates but full distributions for each parameter. Because we modelled an astrophysical phenomenon that in reality must have a single solution (and quantum orbits are not fashionable), we tried to identify a unique orbital configuration by treating the model uncertainty and handling it explicitly. Xiao et al. (2024) assessed that for astrometry data, two data points (i.e., HIPPARCOS and GDR3) are insufficient to determine the inclination, resulting in a bi-modal solution. A third point (from GDR2) resolved this conundrum for their orbital fitting. We diagnosed this phenomenon as the origin of our multimodality. For this particular system, GDR2 leaves little information (poor fit; see Fig. 10), leaving us with a bimodal solution.

Consequently, we tried to fit the modes separately and to quantify the odds. Instead of slicing the parameter space arbitrarily, we set a small normal prior on each possible inclination, with 2K₁~N(0.92,0.3), and 2K₂~N(2.29,0.3). As a result, $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=0.15$ marginally favoured 2K₂ (see the parameter summary for both models in Table D.4). We took the 1K model as our solution, and we further discuss this system in Sect. 5. We also took note of the third and second highest RV periodogram peaks at 403.8 d (which might correspond to the yearly alias with P₁) and 57.7 d, which subsequently disappear in the residuals of 1K and 2K.

Fig. 9 HIP 39330 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.

4.5 HIP 98599

For HIP98599 (or HD 189627), the GLS (see Fig. A.3) reveals a single peak dominating the RVs at 2695 d, and the FWHM GLS has highly significant peaks surrounding it, at 1819.0 d and 3953.9 d. The structure is similar to that of the FWHM in the BIS, but without significant power. The correlogram shows moderate to strong correlations (|ρ| ∈ [0.4,0.8[) between the RVs with FWHM and BIS: ρ_RV,FW=0.64, ρ_RV,BIS=0.51. Similarly to HIP 8923, both correlations might be artefacts due to the uncertainty differences between datasets: the correlations between the RV uncertainty and indices are strong to very strong ρ_RVe,FW=0.96, and ρ_RVe,BIS=0.62. This is further supported by the correlations in the HARPS datasets alone,ρ_RVHAR,FW=-0.26, and ρ_RVHAR,BIS=-0.11, which shift from strongly positive to weakly negative.

The EMPEROR RV-only analysis settles at 1K with $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-144.09\pm 0.14$, with $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=47.50$ over the 0K, with $P_1={2580}_{-90}^{+80}$ d. The RV+astrometry analysis settles at 1K as well, with $P_1={2656}_{-16}^{+40}$ d (see Figs. 11 and 12), with $\ln \hat{\mathcal{Z}}{\mid }_{1\mathcal{K}}=-543.01\pm 0.20$, and $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=61.12$ over the 0K model, a BF difference of 13.62 with the RV-only model. This result also matches the main GLS peak in the RVs. Furthermore, examining the residuals of the GLS yields no periods above the 99.9% FAP line (see Fig. A.3), the most prominent ones (≥ 99.0%) at 7.2 d, 25.5 d and 155.9 d. The residual correlogram also flattens, with ρ_RV,FW = –0.09, ρ_RV,BIS=0.17.

Fig. 10 HIP 39330 best fit for the 2K1 model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

5 Discussion

The EMPEROR framework, with the addition of astrometric differencing modelling (Feng et al. 2019, 2021) introduced in Sect. 2.1, proves to be a useful tool for identifying and characterising cold giants. In our analysis of this first batch from the extended CHEPS (totalling an RV baseline of ∼16 yr), we validated our tool by refining orbital parameter estimates for HIP 21850 and fully characterised five new exoplanets: a warm Jupiter (HIP 10090c), four Jupiter analogues (HIP 8923b, HIP 10090b, HIP98599b, and HIP39330b), and a warm-Saturn candidate (HIP39330c). They are summarised in Table 2. In this section, we quantify how the inclusion of astrometry data boosts our confidence in the planetary solutions via Bayesian evidence and helps us to derive tighter uncertainties for the orbital parameters. Then, we discuss the multi-modality problem presented by HIP 39330, and we address how including small datasets affects our inference. Finally, we place the CHEPS planets of this study in context with the overall exoplanet population, and we highlight future opportunities.

Fig. 11 HIP98599 RVs phase-folded at P1 = 2656-40 d as circles coloured per instrument with the 1K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.

5.1 Improved model evidence

Adding astrometry resulted in periods that were consistent with those found by RVs alone, and it also improved the Bayes factor between models. Defining this increase as $\mathrm{\Delta }\mathrm{B}\mathrm{F}=\mathrm{\Delta }\ln \hat{\mathcal{Z}}{\mid }_{\mathrm{R}\mathrm{V}+\mathrm{A}\mathrm{M}}-\mathrm{\Delta }\ln \hat{\mathcal{Z}}{\mid }_{\mathrm{R}\mathrm{V}}$, we found a wide range of model comparison improvements. In HIP 8923, we have ∆BF = 58.58, whereas in HIP 21850, we have ∆BF = 1.70. To further analyse this increase, we briefly introduce two metrics for our data, a phase-coverage metric Φ_c(P) ∈ [0,1], and a baseline-coverage metric β_c ∈ R⁺. For a given test period P, we mapped to phase each observation time t_i, with i = (0,...,n) for n observations as

$${\varphi }_i\equiv \left(\frac{t_i-t_0}{P}\right)(mod1).$$(7)

Then, we sorted the phases and computed the circular gap g_i,

$$g_i=\{\begin{array}{l}{\varphi }_{i+1}-{\varphi }_i \mathrm{for} \mathrm{i}=1,\dots ,\mathrm{n}-1\\ 1-{\varphi }_n+{\varphi }_1 \text{for i=n}\end{array}.$$(8)

Each g_i reveals the observational gap in phase. With Φ_c ≡ max(g_i), we have a metric for phase incompleteness, where Φ_c → 0 translates into a full period coverage (impossible for a finite sample), and Φ_c → 1 to no phase coverage. We computed Φ_c for a period grid [0, 20 000] for each data type and system, obtaining a periodogram-like figure (see Fig. C.2), where peaks mark the periods that have the least coverage. Additionally, with the baseline coverage B(t) (how many orbital cycles the baseline spans), we calculated the baseline-coverage improvement ratio B_c as

$$B(t)=\frac{t_{\mathrm{m}\mathrm{a}\mathrm{x}}-t_{\mathrm{m}\mathrm{i}\mathrm{n}}}{P},B_c=\frac{B(t_{\mathrm{R}\mathrm{V}+\mathrm{A}\mathrm{M}})}{B(t_{\mathrm{R}\mathrm{V}})}.$$(9)

We used these metrics because they capture the fraction of the orbital phase that is spanned by observations (Φ_c) and the extent to which astrometry extends the temporal baseline (B_c) in an intuitive way.

For HIP 8923, we have ∆BF = 58.58, a tremendous information jump by the inclusion of the astrometry data. This is easily explained: for our period of interest (P₁ ≈ 5162 d), Φ_c decreases considerably with the addition of astrometry data, from Φ_c,RV = 0.468 to Φ_c,RV+AM = 0. 194. This can be seen by comparing the height (Φ_c) of the blue line (RV) against the purple line (RV+AM), across the green vertical bar (P₁ ≈ 5162 d) in Fig. C.2). A summary of comparative values can be found in Table 4. The evidence uncertainties are all ≈ 0.2, and the ∆BF uncertainty is ≈ 0.3. The exact evidence uncertainties can be found in Sect. 4.

At the other extreme, HIP 21850 has the smallest ∆BF = 1.70, which is well explained by Φ_c barely increasing with astrometry ∆Φ_c = 0.016 (in Fig. C.2, our period lies in a peak), and the lowest baseline-coverage ratio B_c = 1.687.

HIP 10090 has the second largest increase ∆BF = 16.30. While not directly explained by coverage metrics, a secondary peak appears in the RV posteriors at 395 d, with ∆likelihood ∼6 from the maximum, which is dispelled with astrometry. In contrast, while it has good coverage metrics, the addition of astrometry creates a multimodal solution for HIP 39330, penalising the overall model likelihood, with the second lowest ∆BF = 3.04. The phase coverage provided by the astrometry in Fig. C.2 shows that our long period resides right in the middle of a red peak Φ_c,AM = 0.723, meaning that astrometry does not provide good coverage at all. The short period (vertical orange line), while being almost completely covered Φ_c,AM ~ 0, is almost imperceptible in the astrometry data (see the small wobble in Fig. 10). It is worth mentioning that the inclusion of astrometry, while mostly uninformative for the 2K model, does indeed help with the first signal. When we compare the BF increase between 0K and 1K, the increase is ∆BF = 9.30.

Fig. 12 HIP 98599 best fit for the 1K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

5.2 Improved planet parameter estimates

Most of our systems present improved parameter estimates. For example, the mass determination helped us to secure HIP 10090c from an ambiguous Saturn m sin i = ${0.39}_{-0.03}^{+0.02}{M}_{\mathrm{J}}$ to a bona fide Jupiter-mass planet M = ${0.85}_{-0.12}^{+0.03}$ M_J. For HIP 21850, we have tighter period constraints $P_1={2534}_{-11}^{+22}\to {2539}_{-3}^{+1}$ d, and for the second period $P_2=15{900}_{-4800}^{+2500}\to 11{320}_{-940}^{+490}$ d, the uncertainties are reduced by almost an order of magnitude, as are the near-coplanar orbits ∆I = 7.2° for this system. For HIP 98599, the period uncertainty range is reduced by a third $P_1={2580}_{-90}^{+80}\to {2656}_{-16}^{+40}$ d.

For HIP 39330, the situation is awkward and effectively leaves us with two astrometric solutions. This means that our data are genuinely insufficient for identifying a unique orbital configuration (see Fig. C.1), bifurcating the second orbital period into two distinct solutions, P_2,1 ≈ 255 d and P_2,2 ≈ 260 d, and leaving us unable to resolve the inclination degeneracy. We fitted each mode separately with normal priors on the inclination I_1,1 ~ N(0.92,0.3), and I_1,2~N(2.29,0.3) and quantified their odds, which are virtually equal. While there is enough evidence to confidently claim a second signal, we cannot confidently characterise it. Conservatively, we prefer adopting the simpler 1K model, until further astrometric data are available to break this degeneracy. The parameter summaries for each mode can be found in Table D.4.

5.3 Sufficient RV data

Our need to split RV datasets into subsets as small as N = 3 (see Table C.1) raises the question whether the information obtained by so few points is worth the additional offset parameter. We examined HIP 10090 in detail, which is composed of a grand total of 57 RVs-12 COR07, three COR14T, six COR14, 32 H, and four H15 (datasets defined in Sect. 3). Out of these datasets, COR07, COR14T, and H15 do not provide additional baseline. Furthermore, their RV phase coverage is shared with other datasets. We compared the full RV phase incompleteness against the RVs without each of these sets (which we denote with the sub-indices −C07, -C14T, and−H15), obtaining for the long period $P_1={2960}_{-100}^{+120}$ the values Φ_c,RV = 0.163, Φ_c,-C07 = 0.213, Φ_c,-C14T = 0.163, and Φ_c,-H15 = 0.163, and for the short period Φ_c=0.136, 0.158, 0.136, and 0.141 for all combined datasets, and -COR07, -COR14T, and -H15, respectively. With the exception of COR07 for the long period, none of these sets has a unique phase coverage, which is reflected in the phase-coverage metric, which remains at 0.163 when datasets are removed. Several interesting points arise when we examine the evidence (see Table 5). Out of these subsets, COR07 drives the most evidence, presenting the highest evidence gap with respect to the full RVs ∼45 for 1K and 2K. As expected, the solution loosens for the long period P₁ = 2920-120 d, while remaining within the uncertainties below unity for $P_2={321.38}_{-0.23}^{+0.69}$ d. Nonetheless, this information does not contribute to the model comparison, with similar $\mathrm{\Delta }\ln \hat{\mathcal{Z}}$ to the full RVs $\mathrm{\Delta }\ln \hat{\mathcal{Z}}=27.22\approx 27.65$. A similar phenomenon occurs for H15, but is driven by the high precision of the data instead of their quantity, where the model comparison stays similar ($\mathrm{\Delta }\ln \hat{\mathcal{Z}}=27.22$), with the long-period uncertainty increasing when removing H15 $P_1={3250}_{-230}^{+130}$ d, and $P_2={320.71}_{-0.04}^{+1.04}$ d staying similarly constrained. On the other hand, and perhaps unexpectedly, the three measurements by COR14T, while contributing the least to the overall evidence of 2K with an almost marginal difference of 1.43 against the full RV, have a crucial role in the model comparison by decreasing the $\mathrm{\Delta }\ln \hat{\mathcal{Z}}$ to 11.88 when removed (from 27.22). These three modest measurements have the peculiarity of overlapping with GOST data, driving the inclination and the astrometric offset posteriors to a tighter solution.

5.4 Cold Jupiters and next steps

Our CHEPS analysis places four systems firmly in the Jupiter-analogue regime (masses of 0.3-13 M_J on orbits of 3-7 AU) and provides a dynamical characterisation for each with modest eccentricities e ≲ 0.19 (see Table 2). Two of the planets, HIP 8923b and HIP39330b, reside very close to Jupiter’s orbital scale (∼5 AU), directly probing Solar System-like architectures. Furthermore, although we were unable to precisely constrain the inclination with our current data, HIP 39330c resides at ∼1 AU, making it a Jupiter-Earth-like system.

Of the more than 6000 exoplanets confirmed to date, only about one-third have true mass estimates (2223 according to Christiansen et al. 2025). Most of these come from transits (1567), so they predominantly sample short to intermediate orbital periods. Figure 13 situates these planets within the broader exoplanet census in the mass-distance plane. The CHEPS Jupiter analogues occupy the transitional swath in which the RV sensitivity declines with increasing a and direct imaging has focused beyond several AU. By combining decades of precise RV monitoring with HIPPARCOS-to-Gaia astrometry, we fill this demographic gap now with dynamical masses instead of minimum masses. This expands the well-vetted sample of true Jupiter analogues, improving constraints for population-level inferences and forward models of direct-imaging yields. Furthermore, the other two planets (HIP 10090c and HIP39330c) occupy a much more sparsely populated region around ∼ 1 AU.

Beyond the characterisation of individual systems, joint RV+astrometry inference is rapidly becoming central to cataloguelevel reliability. In practice, Gaia will turn many long-period RV detections into directly testable hypotheses against Gaia astrometry, either reinforcing the planetary interpretation or exposing a difference that points to additional companions, alternative solutions, or underestimated systematics. This shift will also deliver a step change in mass completeness, moving a large fraction of long-period companions from minimum to true masses and thereby enabling cleaner demographics of cold giants and a more robust separation of planets from brown dwarfs and low-mass stellar companions. Our HIPPARCOS-Gaia differencing and joint-modelling approach already mirrors, in compressed form, the type of information Gaia DR4/DR5 will provide at scale through intermediate astrometry and non-single-star solutions, so the same modelling philosophy and vetting logic will carry over naturally. At the same time, Gaia-based candidate lists will inevitably include astrophysical contaminants (e.g. unresolved binaries), strengthening the case for joint, physically consistent modelling rather than naive catalogue ingestion.

These systems are ripe for targeted follow-up in the near future. Extended RV baselines will continue to refine long-period phase coverage, while upcoming astrometric surveys and data releases such as the Rubin Observatory Data Release 1 on May 2026 (LSST; Hambleton et al. 2023), Gaia ’s DR4 in December 2026 (Brown 2024), the Roman Space Telescope in late 2026-2027 (Holler et al. 2025), or CHES (Ji et al. 2024) will sharpen accelerations and proper-motion anomalies, or even determine orbits on their own with intermediate data. The same RV+astrometry framework that proved decisive here is directly aligned with near-term and mid-term roadmaps for cold giant demographics and precursor target vetting for future missions. Therefore, we encourage the early adoption of available astrometric data in RV characterisation, providing full orbital parameters, as well as enhanced confidence in the modelled solutions, even when RVs alone might appear sufficient. In this context, CHEPS’ metal-rich, quiet FGK targets remain an efficient hunting ground for Solar System analogues that can anchor comparative planet studies across techniques. Building on CHEPS’ 16-year survey of ∼240 stars, already yielding cold giants, we are extending the observational baseline with continued monitoring. This will expand the sample of Jupiter analogues with true masses and full orbits, enabling a robust population-level census of Jupiter analogues.

Fig. 13 True mass vs. semi-major axis of the confirmed exoplanet population from the NASA Exoplanet Archive (Christiansen et al. 2025), coloured by discovery method. The approximate sensitivity curve for GDR4 is shown as a dashed red line, assuming a Sun-like star at 20 pc and a 3σ detection criterion. CHEPS targets from this study are shown as solid black squares. Solar System planets are shown as icons.

6 Conclusions

We extended CHEPS by combining long-baseline RV data with absolute astrometry from the HIPPARCOS and Gaia missions to search for cold, wide-orbiting gas giants. Using the renewed EMPEROR framework enhanced with astrometric differencing, we have analysed five metal-rich FGK hosts and performed Bayesian model comparisons to identify planetary signals and quantify the improvement brought by astrometry.

The analysis confirmed the Jupiter-mass companions for HIP 21850, refined their orbital parameters with CHEPS data by extending the RV-baseline from −16.4 yr to −24.6 yr, and discovered four new Jupiter analogues: HIP 8923b, HIP 10090b, HIP39330b, and HIP98599b. We also identified a new warm Jupiter, HIP 10090c. True masses were derived through the astrometric inclination constraints, yielding estimates from −1 to −10 M_J on orbits between −1 and −7 AU for the new planets (see Table 2). In one system (HIP 39330), the astrometric data led to a multi-modal solution for the second planet; adopting a simpler one planet model is prudent until additional data become available.

Incorporating astrometric data increased the BFs ranging from −1.7 to −58.6, depending on phase coverage, baseline extension, and posterior multi-modality (see Table 4). The custom phase-coverage metrics we introduced showed that astrometry significantly reduces gaps in the orbital phase for cold giants and more than doubles the temporal baseline for some systems. This led to substantial reductions in the period and mass uncertainties, with an improvement of up to an order of magnitude for certain parameters. Even sparse astrometric measurements, when combined with RVs, can break the sin I degeneracy and convert minimum masses into true masses.

The sufficiency test on HIP 10090 illustrated that not all RV subsets contribute equally. While some high-precision but phaseredundant datasets have a marginal effect, others with only a few observations can substantially affect the BF if they sample unique orbital phases. This underscores the need for strategically timed RV observations to complement astrometric data. The inability to constrain the inclination of HIP 39330c also highlights the current limitations of GDR2/GDR3 astrometry for certain orbital periods; Gaia DR4, Rubin DR1, and future missions will provide decisive improvements.

The CHEPS discoveries fill a demographic gap in the exoplanet census, occupying the transition region between RV-dominated detections and direct imaging (see Fig. 13). Establishing a sample of Jupiter analogues with dynamically measured masses will improve population statistics and refine predictions for direct-imaging yields. Continued RV monitoring and forthcoming astrometric releases will tighten the orbital constraints and possibly reveal additional companions. The RV+astrometry framework developed in this work is directly relevant to the nearby targets accessible to ELT coronagraphy and Roman Space Telescope; astrometry provides true masses, orbital geometry, and improved ephemerides that enable scheduling of direct detections and atmospheric characterisation of these cold worlds, which reside at the upper boundary of the planetary mass space.

By combining long-term RV monitoring with absolute astrometry, we demonstrated a powerful avenue for uncovering and characterising cold giant planets on Solar System-like orbits. The framework presented here positions CHEPS and similar surveys to exploit future data releases, paving the way towards a more complete understanding of planetary architectures around nearby stars.

Data availability

Tables with RV measurements with CORALIE per target are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A213. HARPS RVs (Perdelwitz et al. 2024) are available at https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/683/A125. UCLES/AAT RVs can be found in Wittenmyer et al. (2017a). HIRES RVs can be found in Butler et al. (2017). All RV and astrometry products used in this paper can be found in an EMPEROR-friendly format at https://github.com/ReddTea/aa57940-25-data.

Acknowledgements

All benchmarks were performed on a computer with an AMD Ryzen Threadripper 3990X 64-Core Processor with 128Gb of DDR4 3200MHz RAM, limiting the use to 24 physical cores and 24 threads. Typesetting was carried out in Overleaf (https://www.overleaf.com). PAPR and JSJ gratefully acknowledge support by FONDECYT grant 1240738 and from the ANID BASAL project FB210003. We also greatly appreciate the support from the CASSACA China-Chile Joint Research Fund through grant CCJRF2205. Part of this work was supported by the United States Fulbright Fellowship in partnership with the Fulbright Chile Commission, and ANID/Fondo ALMA 2024/No.31240039 “On the origin of Warm-Giant Planets". We thank the anonymous referee for insightful suggestions that enhanced the quality of this manuscript.

References

Appendix A Periodograms

Fig. A.1 HIP 8923 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green vertical region shows P1 = ${5160}_{-240}^{+150}$ d. On the right side there is a table summary of the five highest powers.

Fig. A.2 HIP 10090 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green and orange vertical regions show P1 = ${2960}_{-100}^{+120}$ d and P2 = ${321.82}_{-0.56}^{+0.32}$ d, respectively. On the right side, a table summary with the highest powers.

Fig. A.3 HIP98599 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green vertical region shows P1 = 2656${2656}_{-16}^{+40}$ d. On the right side there is a table summary of the five highest powers.

Fig. A.4 HIP39330periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green and orange vertical regions show P1 = ${4596.9}_{-154.0}^{+1.7}$ d and P2 = ${257.1}_{-1.4}^{+1.6}$ d, respectively. On the right side, a table summary with the highest powers.

Appendix B Correlograms

Fig. B.1 HIP 8923 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.

Fig. B.2 HIP 10090 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.

Fig. B.3 HIP 98599 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.

Fig. B.4 HIP 39330 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.

Appendix C Discussion

Fig. C.1 HIP 39330 corner plot of selected parameters, displaying a bimodal exploration.

Fig. C.2 Phase-coverage per system per data type. The y-axis shows phase-incompleteness Φc (P) as a function of period, with 0 meaning no gaps, increasing up to 1 with gap-size. Blue solid line corresponds to the RV data phase-coverage, red line to astrometry data, and purple to both combined. Green and orange vertical bars display the best-fit orbital periods.

Appendix D Additional tables

All Tables

All Figures

Fig. 1 HIP 21850 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and the window function. The FAP lines are included for 0.1, 1, and 10% and are plotted as dashed red, dotted purple, and dotted blue lines, respectively. The circles show the periods with the greatest power, coloured by FAP region. The coloured regions correspond to $P_1={2539}_{-3.1}^{+0.8}$ d (green) and $P_2={11320}_{-940}^{+490}$ (orange).

In the text

	Fig. 2 HIP 21850 correlogram. Each panel displays the pairwise relation between RVs, FWHM, CCF BIS, and RVe. The Pearson correlation coefficient is denoted by ρ, and the linear trend is shown as a black line.
In the text

	Fig. 3 HIP 21850 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.
In the text

Fig. 4 HIP 21850 best-fit astrometry for the 2K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

In the text

	Fig. 5 HIP 8923 RVs phase-folded at $P={5160}_{-240}^{+150}$ d as circles coloured per instrument with the 1K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right of each plot, we show RV distribution histograms.
In the text

Fig. 6 HIP 8923 best fit for the 1K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

In the text

	Fig. 7 HIP 10090 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right of each plot, we show RV distribution histograms.
In the text

Fig. 8 HIP 10090 best fit for the 2K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

In the text

	Fig. 9 HIP 39330 RVs as circles coloured per instrument with the 2K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.
In the text

Fig. 10 HIP 39330 best fit for the 2K1 model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

In the text

	Fig. 11 HIP98599 RVs phase-folded at P1 = 2656-40 d as circles coloured per instrument with the 1K model imposed as a black line. At the bottom, we plot the RV residuals, and in the right corner of each plot, we show RV distribution histograms.
In the text

Fig. 12 HIP 98599 best fit for the 1K model. Panel a: best-fit astrometric orbit, the black cross is the system barycentre, the grey line connects it to the periapsis, and the grey dotted line joins ascending and descending nodes. Reference epochs are shown as diamonds with their position and proper motion uncertainty as hourglass-shaped shaded regions. Green to blue coloured rings show GOST data projected with the best-fit. Solid circles and slopes show the best-fit position and proper motion induced by the companion. Panel b: zoom-in for the Gaia and GOST model (grey rectangle in panel a). Panel c: residuals for the HIPPARCOS abscissa as solid grey circles, where multiple measurements per epoch have been binned and coloured by time.

In the text

Fig. 13 True mass vs. semi-major axis of the confirmed exoplanet population from the NASA Exoplanet Archive (Christiansen et al. 2025), coloured by discovery method. The approximate sensitivity curve for GDR4 is shown as a dashed red line, assuming a Sun-like star at 20 pc and a 3σ detection criterion. CHEPS targets from this study are shown as solid black squares. Solar System planets are shown as icons.

In the text

Fig. A.1 HIP 8923 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green vertical region shows P1 = ${5160}_{-240}^{+150}$ d. On the right side there is a table summary of the five highest powers.

In the text

Fig. A.2 HIP 10090 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green and orange vertical regions show P1 = ${2960}_{-100}^{+120}$ d and P2 = ${321.82}_{-0.56}^{+0.32}$ d, respectively. On the right side, a table summary with the highest powers.

In the text

Fig. A.3 HIP98599 periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green vertical region shows P1 = 2656${2656}_{-16}^{+40}$ d. On the right side there is a table summary of the five highest powers.

In the text

Fig. A.4 HIP39330periodograms. Shown from top to bottom: RVs, model residuals, FWHM, BIS, and window function. FAP lines for 0.1, 1, and 10%, in dashed red, dotted purple, and dotted blue, respectively. Circle markers show the periods with greatest power, coloured by FAP region. Green and orange vertical regions show P1 = ${4596.9}_{-154.0}^{+1.7}$ d and P2 = ${257.1}_{-1.4}^{+1.6}$ d, respectively. On the right side, a table summary with the highest powers.

In the text

	Fig. B.1 HIP 8923 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.
In the text

	Fig. B.2 HIP 10090 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.
In the text

	Fig. B.3 HIP 98599 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.
In the text

	Fig. B.4 HIP 39330 correlogram. Each panel displays the pairwise relation between RVs, full-width half-maximum (FWHM), CCF bisector inverse slope (BIS), and formal RV uncertainty (RVe). Pearson correlation coefficient is denoted by ρ, with the linear trend as a black line.
In the text

	Fig. C.1 HIP 39330 corner plot of selected parameters, displaying a bimodal exploration.
In the text

	Fig. C.2 Phase-coverage per system per data type. The y-axis shows phase-incompleteness Φc (P) as a function of period, with 0 meaning no gaps, increasing up to 1 with gap-size. Blue solid line corresponds to the RV data phase-coverage, red line to astrometry data, and purple to both combined. Green and orange vertical bars display the best-fit orbital periods.
In the text