© The Authors 2026

1 Introduction

In the study of large sub-stellar companions, one of the main questions is whether the object can be classified as an extrasolar planet or as a brown dwarf (BD). In the transition regime, roughly between 10 and 20 M_J, providing a definitive answer to this question can be very challenging. Traditionally, the threshold between the two types of objects has been set at 13 M_J, which is considered to be the limiting mass for the fusion of deuterium (Boss et al. 2003). However, more recent studies have highlighted that the lower-mass limit for deuterium burning might vary between 11 and 16 M_J, depending on the object composition and other factors (e.g. Mollière & Mordasini 2012; Schneider 2018). Other studies have suggested using the formation mechanisms to distinguish the two populations of companions (e.g. Chabrier et al. 2014), defining the companions formed via gravitational instability as BDs (e.g. Boss & Kanodia 2023), while exoplanets are considered to be those formed via core-nucleated or pebble-assisted accretion within a disc (e.g. Mordasini et al. 2012b). Unfortunately, these definitions are difficult to apply with respect to observed companions because the formation history is difficult to infer from the observed properties, since different formation models can produce similar outcomes. Moreover, the number of currently known sub-stellar companions on the transition regime between BD and GP with precisely measured orbital parameters is still very limited. Most of the currently known long-period objects come from radial velocity (RV) surveys, which means that the true mass is unknown, making it extremely difficult to unveil their true nature.

The advent of ESA’s Gaia mission (Gaia Collaboration 2016) provides a promising new avenue for overcoming these limitations. Gaia’s astrometric capabilities enable precise measurements of stellar positions and motions, making it a powerful tool for detecting the gravitational influence of sub-stellar companions, especially in intermediate to long-period orbits. In particular, astrometric data from Gaia complement RV measurements by offering an independent method for determining the true mass and orbital parameters of companions. The extended Gaia mission now holds the potential to detect GPs and low- mass BDs out to distances of 3–5 AU, for stars of 0.3–1 M_⊙, a critical regime for understanding the transition between these two classes of objects. The first astrometric results from Gaia (Gaia Collaboration 2023b), including detections of sub-stellar companions, are already providing valuable insights into the demographics of these systems, with a few thousand orbits of sub-stellar companion candidates having been published as part of Gaia Data Release (DR) 3 (Gaia Collaboration 2023a). This unprecedented dataset, combined with complementary RV data, sets the stage for a new phase in exoplanet research, offering a more complete view of system formation and evolution.

The first two confirmed Gaia exoplanets were transiting hot-Jupiters detected in the photometry, Gaia-1b and Gaia- 2b (Panahi et al. 2022). Subsequently, Sozzetti et al. (2023) announced the confirmation of the first astrometric Gaia candidate, Gaia-3b, a super-Jupiter on a 300 d period orbit. However, later internal investigation of the Gaia pipeline revealed that the time-series for this target suffered from a software bug and that the orbital solution published in Gaia DR3 was a false positive¹. Since the RV planetary signal has been confirmed, Gaia-3b still represents a Gaia-enabled discovery, as the target would likely not have been monitored without its initial identification dr3-known-issues#FalsePositive as an astrometric candidate². Finally, Stefánsson et al. (2025) announced the confirmation of two more astrometric Gaia planet candidates, Gaia-4b and Gaia-5b, identified as a massive planet (11.8 M_J) at P = 571.3 d and a low-mass BD (20.9 M_J) at P = 358.58 d. These companions were confirmed by combining RV measurements from the Habitable-zone Planet Finder (HPF), NEID, and FIES spectrographs.

In this context, we present results from an intensive high- precision RV monitoring campaign of another Gaia astrometric planet candidate, Gaia-ASOI-009 (aka HD 128717). The object is part of the publicly available Gaia candidate exoplanet list³. The Gaia DR3 archive reports a value of the renormalised unit weight error (RUWE)⁴ statistic of 1.32, formally below the threshold of 1.4 adopted for systematic attempts to perform orbital fits to Gaia DR3 astrometric time-series. The Gaia DR3 solution type, ‘OrbitalTargetedSearch’, corresponds to one of the orbital solutions obtained as part of an analysis performed on a supplementary external input list, as described in detail in Holl et al. (2023) and Gaia Collaboration (2023a). We observed HD 128717 within the context of the Global Architecture of Planetary Systems collaboration (GAPS, Covino et al. 2013; Desidera et al. 2013) with the HARPS-N spectrograph (Cosentino et al. 2012) at Telescopio Nazionale Galileo (TNG). The target was included in the target list of the original GAPS programme, with a tailored observing strategy designed upon release of the Gaia DR3 orbital solution. HD 128717 was then included in the ongoing large programme: The Great HARPS-N hunt for superEarths and Neptunes interior to outer giant planets detected by Gaia. This programme is focussed on confirming astrometric long-period exoplanets and BDs and, through subsequent high- cadence follow-up, as well as on searching for inner low-mass planets in the observed systems. A full description of the programme will be provided in an upcoming paper (Barbato et al., in prep.). Here, we present a detailed analysis of the RV dataset collected within the programme, as well as a combined analysis with the available Gaia DR3 and Hipparcos astrometry to confirm the orbital solution, along with a set of simulations to confirm the planetary nature of the candidate.

In Sects. 2 and 3, we present the new stellar parameters derived for the target from the HARPS-N spectra and the RV data collected for this study, respectively. In Sect. 4, we present the Keplerian fits of the RV time series and the combined analysis of RVs and astrometric data. In Sect. 5, we present our simulations of Gaia observations. In Sect. 6, we combine the PMa astrometric information with our RV data to perform a joint modelling of the orbit. In Sect. 7, we discuss our results and the follow-up imaging observations we performed. Finally, in Sect. 8 we summarise our findings.

2 Spectroscopic HARPS-N data

As part of the GAPS RV programme, HD 128717 was observed with HARPS-N from BJD = 2459735.6 (5 June 2022) to BJD = 2460686.7 (11 January 2025). In total, 106 spectra were collected over a time span of 951 days. These spectra were reduced with version 3.7.1 of the HARPS-N DRS (data reduction software) pipeline, maintained by the Italian Centre for Astronomical Archive (IA2)⁵. All observations were collected with an integration time of 600 s, to average out potential short-term periodic oscillations of the star (Dumusque et al. 2011) and to obtain a sufficient per-pixel signal-to-noise ratio (average S/N = 40 at 550 nm).

The RV measurements were derived using the Template Enhanced Radial velocity Reanalysis Application (TERRA, Anglada-Escudé & Butler 2012). The resulting RV time series has a median error of σ_RV = 1.6 m s⁻¹ and a rms of 196.7 m s⁻¹. The RV dataset is shown in Fig. 1. The complete time series of all spectroscopic data used in this work is available at the CDS. This figure also displays the expected RV signal corresponding to the best-fit orbital solution for HD 128717 b from Gaia DR3 (Gaia Collaboration 2023a). The DR3 orbital parameters are listed in Table 1. The table lists (along with the parameters available in the gaiadr3.nss_two_body_orbits and gaiadr3.binary_masses tables at the Gaia archive)⁶ the derived value of the Campbell elements, along with the values of the minimum mass, m_b sin i, and RV semi-amplitude, K_b.

As we can see, there is a clear inconsistency between the predicted planetary signal and the HARPS-N observations, both in amplitude and phase. It is worth noticing that the orbital period listed in Table 1 corresponds to the 34 months time span of Gaia DR3. This might suggest an underestimated period due to insufficient phase coverage and we discuss this discrepancy in detail in the following sections.

Fig. 1 Upper panel: HARPS-N RV time series of HD 128717 (blue points), compared with the predicted RV signals from the Gaia orbital solution (red dashed line). Lower panel : HARPS-N log R′HK time series.

3 Stellar properties

We derived the effective temperature (T_eff), the surface gravity (log g), the microturbulence velocity (ξ), and the iron abundance ([Fe/H]) using the equivalent width method to the coadded spectrum of the target (see Biazzo et al. 2022 for details). We considered the MOOG code (Sneden 1973; version 2019), the iron line list by Biazzo et al. (2022), and the ATLAS9+ODFNEW (Castelli & Kurucz 2003) grids of model atmospheres. Briefly, T_eff was measured by imposing the excitation equilibrium of Fe I lines, log g through the ionisation equilibrium between Fe I and Fe II lines, and ξ was obtained by removing the trend between the Fe I abundances and the reduced equivalent width REW = log(EW/λ). From this, we found T_eff = 6320±40K, log g = 4.37 ± 0.13, and [Fe/H] = 0.16 ± 0.08.

Moreover, adopting the macroturbulence velocity of 4.7 km s⁻¹ by Doyle et al. (2014), we also derived the projected rotational velocity v sin i_⋆ = 6.1 ± 0.5 km s⁻¹ via the spectral synthesis technique applied to three different wavelength regions of the co-added spectrum (see Biazzo et al. 2022 for the procedure).

We also detected the presence of the lithium line at 6707.8 Å, deriving an equivalent width of 67.0 ± 1.0 mÅ. Adopting our spectroscopic parameters measured above, we obtained a lithium abundance of log(N_Li)_NLTE = 2.78 ± 0.02 dex, after correcting for non-local thermodynamical equilibrium (NLTE) effects (Lind et al. 2009). This high lithium abundance might suggest a young stellar age, discussed further in Appendix C.

We determined the stellar mass, radius, and age by simultaneously modelling the stellar spectral energy distribution (SED) and the MESA Isochrones and Stellar Tracks (Paxton et al. 2015) through the EXOFASTv2 differential evolution Markov chain Monte Carlo tool (MCMC; see Eastman 2017; Eastman et al. 2019 and Naponiello et al. 2025 for more details). To sample the SED, we used the available Tycho-2 B_T and V_T, 2MASS J, H, and K_s, and WISE W1, W2, W3, and W4 magnitudes (see Table 2 and Fig. 2). We imposed Gaussian priors on the T_eff and [Fe/H] from the spectral analysis above, as well as on the Gaia DR3 parallax. From the medians and 16^th – 84^th percentiles of the posterior distributions we found $M_{\star }={1.212}_{-0.068}^{+0.058}$, R_* = 1.248 ± 0.024 R_⊙, and age $t={1.9}_{-1.3}^{+1.9}$ Gyr; the corresponding log g is fully consistent with the^-spectroscopic log g determined from the spectral analysis. All the parameters of the host star are given in Table 2.

4 Spectroscopic data analyses

In this section, we discuss the Keplerian fits of the RV time series of HD 128717. We also present the ancillary analyses of the spectroscopic activity indices, performed to confirm the Keplerian RV signal. We first present a blind analysis of the spectroscopic data, as would be performed to discover and characterise the RV signal from an unknown planetary candidate. We then discuss the combination of the astrometric Gaia DR3 candidate solution and the spectroscopic information for a complete characterisation of the orbital parameters.

The identification of periodic signals in the RV (or activity) time series was performed using a generalised Lomb–Scargle periodogram (GLS, Zechmeister & Kürster 2009). The Keple- rian signal was then fully characterised using a Markov chain Monte Carlo (MCMC) approach, specifically employing the emcee Affine Invariant MCMC Ensemble sampler developed by Foreman-Mackey et al. (2013). We included offset and jitter terms to take into account instrumental effects, and also tested a multi-Keplerian model and a linear RV trend to evaluate the presence of additional companions in the system. The model comparison was performed via the Bayesian information criterion (BIC, Schwarz 1978). The adopted analysis approach is the same as in Pinamonti et al. (2023).

Fig. 2 Stellar spectral energy distribution (SED). The broad-band measurements from the Tycho, 2MASS, and WISE magnitudes are shown in red and the corresponding theoretical values with blue circles. The unaveraged best-fit model is displayed with a black solid line.

Fig. 3 GLS periodograms of the spectroscopic time series. Upper panel : RV. Lower panel : log R′HK.

4.1 Periodogram analyses

Fig. 3 shows the GLS periodograms of the RV and log R′_HK. The RV periodogram identified as the strongest periodic signal at f = 0.000999 ± 0.000025 d⁻¹ (P = 1001 ±25 d), with an amplitude K = 266.1 ± 7.0 m s⁻¹, and a bootstrap false alarm probability, FAP < 10⁻⁴ (Endl et al. 2001). It is worth noticing that, as apparent from Fig. 1, the RV signal present in the time series is highly non-sinusoidal, and thus the sine-wave fitting performed by GLS is a poor match to the data. Nevertheless, the period shows good correspondence with the value from Gaia DR3 (see Table 1).

The amplitude of the signal is much larger than the value derived from the astrometric orbital parameters (i.e. K ≃ 74 m s⁻¹) as observed in Fig. 1. The log R′_HK periodogram shows a significant periodicity at P = 460 d. This however, does not appear to be connected with the RV signal, as the two time series show a very low Pearson (Spearman) correlation, ρ = 0.21 (r_s = 0.25), as shown in Fig. 4. The only other significant peak in the log R′_HK periodogram is at P ≃ 230 d, close to half the first peak period. After removing these two signals, no other significant peak is found.

As an additional test of the influence of stellar activity on the RV time series, we tested the correlation between the log R′_HK and the residuals of the one-Keplerian model (see Sect. 4.2). The correlation appears to be only slightly higher than in the original time series, ρ = 0.26 (r_s = 0.29), as shown in Fig. 5. This again does present evidence for a strong stellar influence on the RV data.

Fig. 4 RV-activity correlation: RV time series as a function of log R′HK.

Fig. 5 RV–activity correlation: RV residuals time series as a function of log R′HK.

4.2 MCMC RV Keplerian fit

We then fit a Keplerian orbit to the HARPS-N RV time series via the MCMC framework described above. We adopted uninformative uniform priors for the Keplerian parameters, semiamplitude, K_b, orbital period, P_b, and time of periastron passage, T_p,b, as listed in Table 3. For the eccentricity and argument of periastron, we followed the $\sqrt{e}\cos \omega$ and $\sqrt{e}\cos \omega$ parametrisation from Eastman et al. (2013). The RV model and residuals are shown in Fig. 6. From this model, we obtained an amplitude of $K_{\text{b}}={321.9}_{-3.4}^{+3.6}{\text{ms}}^{-1}$, a per iod of $P_{\text{b}}={2607}_{-79}^{+84}\text{d}$, and an eccen- +0.0064 tricity of $e_{\text{b}}={0.8061}_{-0.0065}^{+0.0064}$, corresponding to a minimum mass of $M_{\text{b}}\sin i={14.67}_{-0.51}^{+0.50}{M}_{\text{J}}$. The complete posterior distributions of the model are shown in Fig. F.1. These values are significantly different from the Gaia DR3 solution, with the longer period and higher eccentricity resulting in a larger mass, at the threshold between GPs and BDs.

Two outliers can be identified in the RV residuals in Fig. 6, more than 3σ away from the median. These data do not stand out in the log R′_HK time series shown in Fig. 1, so we checked additional indicators to identify the source of these outliers. Checking other CCF activity indicators, the second epoch (BJD=2459844.32) stands out as an outlier in the FWHM and BIS time series (see Appendix A), while the first (BJD=2459830.33) was taken during a night of bad weather and very variable seeing. Thus, we decided to remove both these epochs in all the following analyses to avoid contaminations in the subsequent tests. We then performed a GLS periodogram on the RV residuals to test the presence of additional signal, planetary or stellar, in the time series. As shown in Fig. 7, no significant signal was found.

Fig. 6 Phase-folded best-fit RV one-Keplerian model, corrected for the instrumental offset (top panel) and the RV residuals (bottom panel). The orange dots represents the outliers discussed in Sect. 4.2, while the grey- shaded area represent the phase coverage of Gaia DR3 observations.

Fig. 7 GLS periodogram of the RV residuals of the Keplerian model.

4.3 Additional RV models

Although no significant peak was found in the GLS periodogram of the residuals, we tested a few additional MCMC models to investigate in depth the presence of potential additional signals in the data. We also tested an additional fit to see if the Gaia DR3 orbital solution could be reconciled with the results of our RV Keplerian model.

4.3.1 Long-term trend model

We also tested the presence of additional long-period signals in the data in the form of a linear acceleration term added to the model, $d\cdot (t-\overline{t})$. As shown in Table 3, this model shows a worse BIC than the previous one without the linear trend (ΔBIC = +5). Moreover, the best-fit value for the acceleration term is not significantly different from zero. For these reasons, we decided to discard this model and adopt as final model the previous one.

4.3.2 Two-Keplerian model

As an additional test on the presence of short-period companions in the system, we run an MCMC two-Keplerian model, with wide priors for the second Keplerian period, P_c ∈ 𝒰(1.5,100) d (model not shown). The MCMC did not converge to any significant periodicity, and the two-Keplerian model is statistically disfavoured with respect to the one-Keplerian one (ΔBIC = +5).

5 Gaia simulations

To further study the discrepancy between the RV and astrometric solutions, we performed a suite of numerical simulations. First, we produced a set of realistic Gaia astrometric time series, injecting possible realisations of a Keplerian signal corresponding to the RV orbital fit. We then studied the orbital solutions recovered by fitting these synthetic astrometric data. This process was meant to help us explain why the Gaia DR3 solution was so different from our findings and to further study whether Gaia DR4 could help us confirm our results.

The setup of this simulation is the same as in Sozzetti et al. (2023) and Ruggieri et al. (2024). We used GOST⁷ to obtain the Gaia observation information (observation time, scan angle, and parallax factor) and took the five-parameter stellar motion solution from Gaia DR3. We took the RV orbital parameters from gost/index.jsp our HARPS-N solution, drawing them from a Gaussian distribution. We generated the longitude of the ascending node, Ω, drawn from a uniform distribution, [0, π], and the inclination, i, drawn from a uniform distribution of cos i ∈ [−1, 1].

5.1 DR3 simulations

We simulated the Gaia DR3 time series obtaining from GOST all the Gaia observations of HD 128717 from July 2014 to May 2017 (Gaia Collaboration 2021). As measurement error, we adopted σ = 150 μas, which is representative of the DR3 astrometric uncertainties for a star of magnitude G = 8.2 (Holl et al. 2023).

We produced 100 realisations of the astrometric time series and fit each one with a partly linearised model (Holl et al. 2023) implemented with emcee. The best-fit values of the inclination, i, and astrometric semi-major axis, a₀, compared with the injected values, are shown in Fig. 8. We can see that although there is a partial correspondence between the injected and recovered values of i (lower panel Fig. 8), the semi-major axis, a₀, is almost never recovered correctly (see upper panel Fig. 8). This also implies that the orbital period is very hard to constrain from the DR3 astrometric data. This is highlighted in Fig. 9, where the distributions of period and eccentricity are compared with the values from the RV and Gaia DR3 orbital solutions. We can see that the distributions do not resemble the RV values, although they were used to generate the synthetic time series and are very wide. In particular, the eccentricity peaks at lower values, close to the value of e = 0.39 obtained by Gaia DR3. If we compute how many times the fitted value is within 1σ from the RV value of P and e, we get six and three times out of a total of 100 simulations, respectively. For the Gaia DR3 solution, we fall within 1σ 10 and 42 times, respectively. Thus, we can see that even if there is no clear bias towards the Gaia DR3 orbital solution, it can be easily recovered from the astrometric time series thanks to the very wide distribution of the possible solutions. This is not unexpected, since the RV data clearly indicate an orbital period much longer than the duration of DR3 observations and, thus, the signal is very difficult to constrain due to the poor orbital coverage.

Fig. 8 Astrometric simulation results for DR3. Upper panel: simulated (red) and derived (blue) astrometric semimajor axis vs injected orbital inclination; Lower panel: derived inclination vs the simulated value.

Fig. 9 Distribution of the recovered values of period, P, on the left, and eccentricity, e, on the right, in the DR3 astrometric simulations. The green lines and shaded regions mark the best-fit values and uncertainties from the Gaia DR3 orbital solution, while the orange lines and regions represent our RV solution.

5.2 DR4 simulations

Next, we produced a set of 100 Gaia DR4 time series with the GOST observations from July 2014 to December 2019. We adopted internal uncertainties of 80 μas, which is expected to be a representative value for similar stars in DR4 (Brown 2024). We can see from Fig. 10 that the correspondence between the injected and recovered values is much better with the DR4 time series, with respect to the DR3 simulations discussed previously. Although some uncertainty on the semi-major axis remains, the orbital inclination is recovered with great accuracy. In Fig. 11, we can see how the distribution of the output values of P and e closely peaks to the RV orbital solution used to produce the datasets. These results show how the current DR3 time span is insufficient to correctly recover the Keplerian orbit, while the next data release is expected to have sufficient coverage to accurately recover the companion signal.

Fig. 10 Astrometric simulation results for DR4. Upper panel: simulated (red) and derived (blue) astrometric semimajor axis vs injected orbital inclination; Lower panel: derived inclination vs simulated value.

Fig. 11 Distribution of the recovered values of period, P, and eccentricity, e, in the DR4 astrometric simulations. The green lines and shaded regions mark the best-fit values and uncertainties from the Gaia DR3 orbital solution, while the orange lines and regions represent our RV solution.

6 Combined fit of RVs and absolute astrometry

As reported in cross-calibrated Hipparcos-Gaia catalogues of astrometric accelerations (Brandt 2021; Kervella et al. 2022), HD 128717 is clearly identified as an accelerating star, with a S/N of the proper motion anomaly (PMa) of ~ 9 at the mean epoch of Gaia DR3. We performed a combined PMa+RV model to the Hipparcos-Gaia absolute astrometry of HD 128717 from Brandt (2021) and to the offset-corrected RV dataset in order to determine the possible values of orbital inclination, i, longitude of the ascending node, Ω, and mass ratio, q, which can be directly constrained by the PMa astrometric data. We adopted uninformative priors for the model parameters derived by the RV-only analysis and uniform priors on cos(i), Ω, and q.

The combined analysis utilises a differential evolution MCMC (DE-MCMC) algorithm (Ter Braak 2006; Eastman et al. 2013), with the PMa model constructed by averaging over the actual observing windows of Hipparcos (utilising the exact time of observation of the target; see van Leeuwen 2007) and Gaia DR3 (based on the accurate representation of the actual Gaia transit times from GOST. The DE-MCMC analysis was run with standard prescriptions on the burn-in phase and convergence conditions based on the Gelman-Rubin statistics (e.g. Eastman et al. 2013). The medians of the posterior distributions of the model parameters constituted the nominal best-fit values, with 1σ uncertainties determined evaluating the ±34th percentile intervals of the posteriors (for additional details, see Sozzetti 2023).

Table 4 reports the best-fit orbital solution for HD 128717 B. Usually, applications of the PMa technique (e.g. Li et al. 2021; Sozzetti 2023) return bimodal distributions of the inclination angle, due to the difficulty in separating prograde from retrograde orbits. In this case, instead of being bimodal, the inclination distribution has a single value of $i_{\text{b}}={130.3}_{-1.9}^{+1.6}$ (with a corresponding ${\text{Ω}}_{\text{b}}={131.7}_{-3.4}^{+3.3}\mathrm{deg}$), possibly because the time interval encompassed by the Gaia DR3 observations critically samples the high-velocity phase of the very eccentric orbit of HD 128717 B just after periastron passage (see Fig. 6). We obtained a mass ratio q_b = 0.0148 ± 0.0004 and derive a true mass of M_b = 19.8 ± 0.5 M_J. The combined solution indicates that HD 128717 B’s orbit has a period about 30% longer and a slightly higher eccentricity than that inferred by the RV-only fit. In Figs. 12 and 13, we show the joint posteriors of these parameters and the best-fit orbital solution superposed to the observational data, respectively.

7 Discussion

In light of the aforementioned analyses, we confirmed the Kep- lerian nature of Gaia-ASOI-009 b. It is a sub-stellar companion of mass M_b = 19.8 ± 0.5 M_J (hereafter, Gaia-6 B).

Fig. 12 Joint posterior distributions for inclination, longitude of the ascending node and true mass for HD 128717 B.

7.1 Comparison with Gaia DR3 solution

Thanks to high-precision RV follow-up with the HARPS-N spectrograph of HD 128717, we were able to confirm the presence and sub-stellar nature of the candidate Gaia-ASOI-009 b. Nonetheless, our final Keplerian solution for Gaia-6 B (see Table 4) is significantly different from the one from Gaia DR3 (Table 1). An additional test to see whether the two solutions could be reconciled can be found in Appendix B but, as expected, the evidence from the RV time series strongly points away from the Gaia DR3 solution.

This discrepancy can be explained using the simulations we performed, as discussed in Sect. 5. As an additional test, we produced 100 Gaia DR3 time series using the orbital parameters from the final PMa+RV solution (Table 4) and compared the likelihood of the simulated data given two different set of parameters from the Keplerian model: the set derived from our PMa+RV fit and the one from the Gaia DR3 solution. We find that the two distributions of likelihood are indistinguishable: a Kolmogorov–Smirnov (K–S) test on the likelihood distributions returned a p-value p = 96%. Moreover, we see that in 42/100 cases the likelihood of the DR3 solution is higher than that of our final solution. This proves that the Gaia DR3 time series is not able to constrain the orbit of Gaia-6 B. The reason for this is most probably the orbital period of the sub-stellar companion: the time span of Gaia DR3 observations (i.e. 34 months) is much shorter than the orbital period, P_b = 9.37 yr. Combined with the high eccentricity of the orbit, this results in a strong degeneracy between eccentricity and orbital period, which most likely produced the Gaia DR3 solution.

Fig. 13 Best-fit orbital solution (solid black line) for HD 128717 B, superposed to the calibrated proper motions from Hipparcos and Gaia (top and middle panels) and to the RV time series (bottom panel). The dashed orange lines represent random selections of orbital solutions drawn from the posterior distributions.

7.2 Giant planet-brown dwarf classification

Given its high mass, M_b = 19.8 ± 0.5 M_J, the sub-stellar companion identified in the astrometric and RV data lies above the commonly adopted 13 M_J threshold dividing ‘planets’ and ‘BDs’; however, the exact threshold is somewhat arbitrary (e.g. Stefànsson et al. 2025, and references therein).

In Fig. 14, we show a comparison of Gaia-6 B with the known population of GPs and BDs, as defined in Appendix E. We can see that the mass of the host star is very close to the peak of both the BD and GP distributions. On the other hand, it has a higher-than-average metallicity, close to the peak of the GP distributions. This can lead to some uncertainty on the nature of the sub-stellar companions, depending on the definition adopted to distinguish planets from BDs. Metallicity is commonly seen as a statistical indicator to distinguish between GPs and BDs (e.g. Ma & Ge 2014; Mata Sànchez et al. 2014; Maldonado & Villaver 2017), as it can be correlated with the formation mechanisms: GPs formed via core accretion are more common around metal-rich stars (e.g. Osborn & Bayliss 2020), while BDs formed via gravitational instability are considered to be independent of metallicity (e.g. Maldonado et al. 2019). From this point of view Gaia-6 B might belong to both populations, being located at the joint between the two populations, but its high metallic- ity could suggest formation through core accretion. However, its mass is significantly higher than the expected lower mass limit for deuterium burning (11–16 M_J, Mordasini et al. 2012a), which strongly points towards it being a deuterium-burning BD.

7.3 Search for short-period companions

In Sect. 4.3, we tried to fit additional signals in the RV time series to test the presence of additional objects in the system. We found no significant signal in the periodogram of the RV residuals and all other models showed a higher BIC than the single-Keplerian one; namely, no statistical evidence of a more complex model was found. Nonetheless, small planetary companions could still be undetected in the HARPS-N data. To quantify this, we computed the detection function of the RV data, following the same formalism adopted in Ruggieri et al. (2024) and Sozzetti et al. (2024). The resulting detection map is shown in Fig. 15. It is worth noticing that low-mass super Earths (M_p sin i ≲ 10 M_⊕) would not be detected in our time series, except at the shortest periods. This in spite of the high-precision, high cadence HARPS-N RV time-series, because of the significant levels of RV jitter.

However, due to the high mass and large eccentricity of Gaia-6 B, many orbital configurations would be unstable. This means that we can rule out the presence of planetary companions even below the detectability threshold. We can compute the orbital stability region using Hill’s Criterion (Hamilton & Burns 1992). The region of the period-mass parameter space⁸ in which a planet on a circular orbit would be unstable under Hill’s Criterion is highlighted in orange in Fig. 15. We can conclude that the only possible undetected planetary companions in the inner regions have M_p sin i ≲ 10 M_⊕ and P ≲ 50 d.

7.4 Imaging follow-up observations

The high eccentricity of HD 128717 B can be interpreted as the result of orbital excitation caused by a massive companion at large separation. As the orbital coverage provided by our HARPS-N RV data is too low to robustly detect any long-term trend caused by such a distant companion, we searched for any outer companion to HD 128717 with imaging, a detection technique uniquely suited to discover massive companions on large orbits beyond the reach of RV observations.

We observed HD 128717 with SHARK-NIR (Farinato et al. 2022; Marafatto et al. 2022) and LBTI/LMIRCam (Skrutskie et al. 2010; Leisenring et al. 2012), two infrared imagers mounted on the Large Binocular Telescope (LBT) in Arizona for multi-band simultaneous observations, on April 15 and May 16, 2025. Details of these observations are discussed in Appendix D. In light of the longer duration and better seeing, in the following, we focus on the May 16 observation alone.

No source with a S/N ratio higher than 4 is present on either of the post-processed images, allowing us to conclude that no companions were successfully detected around HD 128717 during our observations. This non-detection can be used to provide robust limits on the mass of any yet-undetected outer companion in the system. Using the AMES-COND models Allard et al. (2003) and the star host age of 1.4 ± 0.3 Gyr (computed as discussed in Appendix C), we derived the mass detection limits for the imaging observations, shown in the bottom panel of Fig. 16. While no low-mass GP companion would have been detectable in our data, our observations would have instead been able to robustly detect any large-separation stellar or BD companion in the system. Specifically, focussing on the mass limit computed for the nominal age of 1.4 Gyr and orbital separations beyond 60 au, LMIRCam observations in the L′ band are able to exclude the presence of any outer companion more massive than ~47 M_Jup, while the SHARK-NIR data in the H band allow for the robust exclusion of any distant companion above ~11 M_Jup.

Fig. 14 Masses of GPs and BDs as a function of stellar host mass and metallicity, for sub-stellar companions with period < 104 d. BDs (M > 13 MJ) are shown as blue points, while GPs (0.1 < M < 13 MJ) are shown as orange points. The red diamonds highlight the position of confirmed Gaia- ASOI sub-stellar companions (and HIP 66074 B Sozzetti et al. 2023). The lower plots show the stellar mass and metallicity distributions of BD and GP.

Fig. 15 Detection map for the HARPS-N RV data of HD 128717. The colour scale represents the detection function. The red point shows the position of HD 128717 B in the parameter space. The orange shaded area is the forbidden region due to the Hill stability criterion.

Fig. 16 Contrast curve and mass limits for the SHARK-NIR (blue) and LMIR-Cam (orange) observations of HD 128717 conducted on the night of May 16, 2025 UT. Top panel: contrast curve. Bottom panel: mass limits, with the thick curve corresponding to the nominal 1.40 Gyr stellar age and the shaded region to the age uncertainty.

8 Conclusion

We have confirmed the sub-stellar nature of the astrometric candidate Gaia-ASOI-009, thanks to a high-cadence high-precision RV follow-up with the HARPS-N spectrograph. In brief, Gaia- 6 B is a high-eccentricity BD, M_B = 19.8 M_J, e_B = 0.85, and it is one of the most eccentric BDs with a precisely measured mass known so far (see Fig. E.1). Moreover, we solved the apparent discrepancy between the orbital measurement from Gaia DR3 and the high-precision HARPS-N RV follow-up. This proves the importance ofRV follow-up of astrometric candidates to confirm and improve the precision on the orbital parameters.

The origin of the high eccentricity of Gaia-6 B remains an unsolved puzzle, since we did not identify any hint of the presence of other companions in the RV, astrometric, and direct imaging data analysed in this study. This is unexpected, since all known high-eccentricity (> 0.9) planetary companions have clearly identified additional companions on wide orbits (Sozzetti et al. 2023, and references therein) and there is a known link between eccentricity and stellar multiplicity (e.g. Moutou et al. 2017). However, the current data are not sufficient for us to conclusively exclude the presence of small BD companions at around a few tens of astronomical units, which might still be sufficient to excite Gaia-6 B eccentricity to the observed level. In the future, additional observations, both additional RVs data to expand the temporal baseline and deeper imaging observations to search for smaller companions, will help tighten the constraints on the presence of undetected companions. Meanwhile, in-depth dynamical simulations might help shed light on the possible mechanisms producing such high eccentricity, which might be more common than expected among long-period GPs and BDs.

Data availability

Full RV and activity HARPS-N data (Sect. 2 and Appendix A) are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/707/A67

Acknowledgements

We acknowledge support from the European Union – NextGenerationEU (PRIN MUR 2022 20229R43BH) and the “Programma di Ricerca Fondamentale INAF 2023”. We acknowledge financial contribution from the INAF Large Grant 2023 “EXODEMO”. MPi acknowledges support from ASI-INAF agreement no. 2025-10-HH.0 “Partecipazione Italiana al Gaia DPAC – Supporto alle attività di responsabilità del team scientifico”. TZi acknowledges support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0, NVIDIA Academic Hardware Grant Program for the use of the Titan V GPU card and the Italian MUR Departments of Excellence grant 2023–2027 “Quantum Frontiers”. L.M. acknowledges financial contribution from PRIN MUR 2022 project 2022J4H55R. We acknowledge the Italian center for Astronomical Archives (IA2, https://www.ia2.inaf.it), part of the Italian National Institute for Astrophysics (INAF), for providing technical assistance, services and supporting activities of the GAPS collaboration. We thank Tom Herbst of MPIA-Heidelberg (D) and the NIRVANA team for sharing with us part of the NIRVANA instrument control SW to control the motorised axis of SHARK-NIR, and NASA and the PI of JWST/NIRCam [26] Marcia Rieke for allowing us to use one of the NIRCam spare detectors as scientific detector for the SHARK-NIR scientific camera. Eventually, we emphasise that observations have benefited from the use of ALTA Center (alta.arcetri.inaf.it) forecasts performed with the Astro- Meso-Nh model. Initialisation data of the ALTA automatic forecast system come from the General Circulation Model (HRES) of the European Centre for Medium Range Weather Forecasts. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program, and of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

Appendix A Additional activity analysis

In addition to the log R′_HK analysis, we studied two other activity indicators, derived from the CCF of the HARPS-N spectra via version 3.0.1 of the HARPS-N Data Reduction Software (DRS, Dumusque et al. 2021): the full width half maximum (FWHM) and bisector span (BIS). The two time series are shown in Fig. A.1. Of the two outliers discussed in Sect. 4.2, one can be clear seen as an outlier in both CCF activity time series.

Fig. A.1 Time series of the HARPS-N CCF activity indicators. Upper panel: FWHM; Lower panel: BIS. The orange dots represents the outliers discussed in Sect. 4.2.

In Fig. A.2, we show the GLS periodograms of the two time series. The FWHM shows a significant peak at P = 422 ± 27 d, close to the one highlighted in the log R′_HK time series (see Fig. 3), while the BIS does not show any significant periodicity above the FAP= 10% level.

Appendix B Gaia DR3 priors model

As an additional test to attempt reconciling the results from Gaia astrometry and HARPS-N RVs, we performed an MCMC fit of the RV time series using Gaussian priors from the Gaia orbital solutions, using the best-fit values and their uncertainties. The adopted values of these priors are listed in Table B.1. It is worth noticing that while some parameters come directly from the orbital solutions published in the Gaia archive, such as those on P_b and T_p,b, the others were derived from the related parameters and corresponding uncertainties (see Sect. 2). One particular case is that of the semi-amplitude K_b, which directly depends on the companion mass, for which no uncertainty was provided in the Gaia DR3 catalogue: for this reason, we did not fix a Gaussian prior on K_b and adopted the uninformative prior from Table 3. The results of this fit are listed in the last columns of Table B.1: there is no significant deviation from the results of the original 1-Keplerian fit with uninformative priors from Table 3. This suggests that the evidence from the RV time series is stronger than the constraint from the astrometric orbital solution, as expected.

Fig. A.2 GLS periodograms of the HARPS-N CCF activity indicators. Upper panel: FWHM; Lower panel: BIS.

Appendix C Stellar age and rotation period

As discussed in Sect. 3, the isochrones estimate of the stellar age points towards an intermediate age, although with large uncertainties, $t_{iso}={1.9}_{-1.3}^{+1.9}$. We consider here additional age diagnostics to further c^-onstrain the system age.

In Sect. 3, we highlight the presence of a significant Lithium line in the spectra of HD 128717. The abundance of Lithium is commonly associated with a young stellar ag, as Lithium depletes over time in main-sequence stars. However, comparing HD 128717 Lithium equivalent width and abundance, and its effective temperature with those of clusters of known age, the parameters of HD 128717 are in good agreement with the distribution of the open cluster NGC 752, which has a measured age of 1.34 Gyr (Agüeros et al. 2018).

Another common proxy for stellar age is the activity level: the mean log R′_HK = −4.947 is quite low. Following activityrotation and gyrochronology relations (Noyes et al. 1984; Mamajek & Hillenbrand 2008), we can derive a rotation period of P_rot = 21 ±4 d and an age of t = 3.7 ±0.4 Gyr. However, this rotation period is incompatible with the measured projected rotational velocity, v sin i = 6.1 ±0.5 km s⁻¹ (in good agreement with the vsin i T_eff relation by Winn et al. (2017)), that combined with the stellar radius R_⋆ = 1.248 ± 0.024 R_⊙ results in a maximum rotation period of P_{rot, max} = 10.3 ± 0.9 d.

To obtain a more precise estimate of the rotation period, we analysed 4 available TESS sectors (14,15,49 and 50) with the Python wrapper juliet⁹ (Espinoza et al. 2019) using a quasi-periodic Gaussian process (GP) kernel implemented with celerite (Foreman-Mackey et al. 2017). The GP modelling converges on a period of 3.9 ± 0.3 d.

If we interpreted this as the stellar rotation period, it would point towards a young age (≃ 300 Myr), in contrast with the other indicators. However, the peak in the GP posterior could correspond to the first harmonic of the stellar rotation, resulting in P_rot = 7.8 ± 0.6 d, which would correspond in a more compatible estimate of gyrochronological age, t_gyro = 900 Myr.

Finally, we were able to obtain another estimate of the age of the system from its kinematic. We can compute the Galactic space velocity components U, V, and W from the proper motion, parallax, and radial velocity of the system. We can get the Gaia DR3 parallax from Table 2 and, from the HARPS-N spectra collected, we can derive the absolute systemic velocity of the star using version 3.0.1 of the HARPS-N Data Reduction Software (DRS, Dumusque et al. 2021): γ_syst = +5.90776 ± 0.00101 km s⁻¹. We do not use the Gaia DR3 proper motion components, as the measurements can be biased by the incorrect modelling of the Keplerian motion. For this reason, we adopt the proper motion from Brandt (2021), which is computed from the difference in position between the Gaia DR3 and Hipparcos catalogues, and therefore should not be affected by the orbital motion¹⁰: μ_ra = −81.992 ± 0.021 mas yr⁻¹, μ_dec = +67.612 ± 0.023 mas yr⁻¹. From these values, following the methodology by Johnson & Soderblom (1987), we obtain: U = −35.876 ± 0.059 km s⁻¹, V = −10.698 ± 0.013 km s⁻¹, W = −3.673 ± 0.005 km s⁻¹. From these Galactic space velocities, we can compute the kinematical age of the system, as defined by Almeida-Fernandes & Rocha-Pinto (2018), as t_kin ≃ 1.4 Gyr. It is worth mentioning that the kinematical age gives only an approximate estimate for singles systems, and it is more informative for statistical samples, and in fact Almeida-Fernandes & Rocha-Pinto (2018) suggest an uncertainty of around 3 Gyr for individual stars. However, it is noteworthy that this value is in good agreement with the previous estimates. Combining all these estimates, we can adopt the average value of the age for HD 128717 t = 1.4 ± 0.3 Gyr.

Appendix D Direct imaging observations and analysis

In both observation epochs SHARK-NIR used its board-band H filter and Gaussian coronagraph with an inner working angle (IWA) of 120 mas, while LMIRCam observations were performed in the L′ band using the vector-apodising phase plate coronagraph (vAPP, Doelman et al. 2021) having an IWA of 246 mas. The coronagraphs IWAs correspond for HD 128717 to orbital separations of ~8 and ~18 au, meaning that HD 128717 B would remain covered by the central occulter even at apoastron and therefore undetectable. The April observations lasted ~40 minutes and where characterised by rapidly variable seeing between 1.09″ and 2.30″, while the May epoch lasted ~1 hour with seeing between 0.90″ and 1.49″.

The SHARK-NIR and LMIRCam data thus collected were reduced using a post-processing procedure based on the angular differential imaging method (ADI, Marois et al. 2006) method and the principal components analysis (PCA Soummer et al. 2012) algorithm. For both instruments, contrast limits were computed by estimating the standard deviation within one-pixel-wide annuli centred on the star, with the resulting contrast curves shown in the upper panel of Fig. 16.

Appendix E BD and GP statistical sample

To compare Gaia-6 B with the known population of BDs and GPs, we collected from the public archives a large sample of well characterised systems with BD and GP companions. First we selected from the NASA Exoplanet Archive¹¹ all objects with mass (or m sin i) > 0.1 M_J and period < 10⁴ d, considering only the default solutions from the archive. Then we selected only the systems with a measured stellar mass and metallicity, resulting in 1138 companions. We then selected from the Encyclopaedia of Exoplanetary Systems¹² all BDs with mass (or m sin i) > 13. M_J, selecting again only companions with period < 10⁴ d and measured mass and metallicity of the stellar host, corresponding to 78 companions. We then merged this sample with the catalogue of BD companions from Gaia DR3 selected by Stevenson et al. (2023) which, after applying the same selection as for the other two catalogues, contains 117 objects. Merging this three samples together, discarding duplicates, we obtain a sample of 1311 companions, with 1118 and 193 GPs and BDs, respectively. Fig. 14 shows the distribution of companion mass as a function of the stellar mass and metallicity, while Fig. E.1 shows the period-eccentricity distribution of GPs and BDs. In both figures are highlighted Gaia-4,-5, and −6 B as well as HIP66074 b.

Appendix F emcee auxiliary plots

Fig. F.1 shows the complete posteriors of the 1-Keplerian model of the HARPS-N RV data, which was adopted as the best RV model.

Fig. E.1 Eccentricity of GPs (orange) and BDs (blue) as a function of the orbital period, for substellar companions with period < 104 d. The red diamonds highlight the positions of the confirmed Gaia-ASOI sub-stellar companions and HIP 66074.

Fig. F.1 Posterior distributions of the fitted parameters of the 1-Keplerian model of the HARPS-N RV time series. The vertical dashed lines denote the median and the 16th and 84th percentiles. The shown parameters are, in order: kb, Pb, Tp,b, eb, ωb, γHARPS-N, and σjit,HARPS-N

All Tables

All Figures

	Fig. 1 Upper panel: HARPS-N RV time series of HD 128717 (blue points), compared with the predicted RV signals from the Gaia orbital solution (red dashed line). Lower panel : HARPS-N log R′HK time series.
In the text

	Fig. 2 Stellar spectral energy distribution (SED). The broad-band measurements from the Tycho, 2MASS, and WISE magnitudes are shown in red and the corresponding theoretical values with blue circles. The unaveraged best-fit model is displayed with a black solid line.
In the text

	Fig. 3 GLS periodograms of the spectroscopic time series. Upper panel : RV. Lower panel : log R′HK.
In the text

	Fig. 4 RV-activity correlation: RV time series as a function of log R′HK.
In the text

	Fig. 5 RV–activity correlation: RV residuals time series as a function of log R′HK.
In the text

	Fig. 6 Phase-folded best-fit RV one-Keplerian model, corrected for the instrumental offset (top panel) and the RV residuals (bottom panel). The orange dots represents the outliers discussed in Sect. 4.2, while the grey- shaded area represent the phase coverage of Gaia DR3 observations.
In the text

	Fig. 7 GLS periodogram of the RV residuals of the Keplerian model.
In the text

	Fig. 8 Astrometric simulation results for DR3. Upper panel: simulated (red) and derived (blue) astrometric semimajor axis vs injected orbital inclination; Lower panel: derived inclination vs the simulated value.
In the text

	Fig. 9 Distribution of the recovered values of period, P, on the left, and eccentricity, e, on the right, in the DR3 astrometric simulations. The green lines and shaded regions mark the best-fit values and uncertainties from the Gaia DR3 orbital solution, while the orange lines and regions represent our RV solution.
In the text

	Fig. 10 Astrometric simulation results for DR4. Upper panel: simulated (red) and derived (blue) astrometric semimajor axis vs injected orbital inclination; Lower panel: derived inclination vs simulated value.
In the text

	Fig. 11 Distribution of the recovered values of period, P, and eccentricity, e, in the DR4 astrometric simulations. The green lines and shaded regions mark the best-fit values and uncertainties from the Gaia DR3 orbital solution, while the orange lines and regions represent our RV solution.
In the text

	Fig. 12 Joint posterior distributions for inclination, longitude of the ascending node and true mass for HD 128717 B.
In the text

	Fig. 13 Best-fit orbital solution (solid black line) for HD 128717 B, superposed to the calibrated proper motions from Hipparcos and Gaia (top and middle panels) and to the RV time series (bottom panel). The dashed orange lines represent random selections of orbital solutions drawn from the posterior distributions.
In the text

Fig. 14 Masses of GPs and BDs as a function of stellar host mass and metallicity, for sub-stellar companions with period < 104 d. BDs (M > 13 MJ) are shown as blue points, while GPs (0.1 < M < 13 MJ) are shown as orange points. The red diamonds highlight the position of confirmed Gaia- ASOI sub-stellar companions (and HIP 66074 B Sozzetti et al. 2023). The lower plots show the stellar mass and metallicity distributions of BD and GP.

In the text

	Fig. 15 Detection map for the HARPS-N RV data of HD 128717. The colour scale represents the detection function. The red point shows the position of HD 128717 B in the parameter space. The orange shaded area is the forbidden region due to the Hill stability criterion.
In the text

	Fig. 16 Contrast curve and mass limits for the SHARK-NIR (blue) and LMIR-Cam (orange) observations of HD 128717 conducted on the night of May 16, 2025 UT. Top panel: contrast curve. Bottom panel: mass limits, with the thick curve corresponding to the nominal 1.40 Gyr stellar age and the shaded region to the age uncertainty.
In the text

	Fig. A.1 Time series of the HARPS-N CCF activity indicators. Upper panel: FWHM; Lower panel: BIS. The orange dots represents the outliers discussed in Sect. 4.2.
In the text

	Fig. A.2 GLS periodograms of the HARPS-N CCF activity indicators. Upper panel: FWHM; Lower panel: BIS.
In the text

	Fig. E.1 Eccentricity of GPs (orange) and BDs (blue) as a function of the orbital period, for substellar companions with period < 104 d. The red diamonds highlight the positions of the confirmed Gaia-ASOI sub-stellar companions and HIP 66074.
In the text

	Fig. F.1 Posterior distributions of the fitted parameters of the 1-Keplerian model of the HARPS-N RV time series. The vertical dashed lines denote the median and the 16th and 84th percentiles. The shown parameters are, in order: kb, Pb, Tp,b, eb, ωb, γHARPS-N, and σjit,HARPS-N
In the text