The GAPS Programme at TNG

The nearby mid-K dwarf HIP66074 was recently identiﬁed as host to a candidate super-Jupiter companion on a ∼ 300day, almost edge-on, orbit, based on Gaia Data Release 3 (DR3) astrometry. Initial attempts at conﬁrming the planetary nature of the signal based on publicly available radial-velocity (RV) observations uncovered an intriguing conundrum: the inferred RV semi-amplitude appears to be a factor of 15 smaller than the one predicted based on the Gaia solution (corresponding to a 7-M Jup companion on a close to edge-on orbit). We present the results of intensive RV monitoring of HIP66074 with the HARPS-N spectrograph. We detected the companion at the Gaia period, but with an extremely eccentric orbit ( e = 0 . 948 ± 0 . 004), a semi-amplitude K = 93 . 9 + 9 . 4 − 7 . 0 ms − 1 , and a minimum mass m b sin i b = 0 . 79 ± 0 . 05 M Jup . We used detailed simulations of Gaia astrometry with the DR3 time-span to show that the conundrum can be fully resolved by taking into account the combination of the initially sub-optimal RV sampling and systematic biases in the Gaia astrometric solution, which include an underestimation of the eccentricity and incorrect identiﬁcation of orbital inclination, which has turned out to correspond to a close to face-on conﬁguration ( i (cid:46) 13 ◦ ). With an estimated mass in the approximate range of 3 − 7 M Jup , we ﬁnd that HIP66074b ( ≡ Gaia-3b) is the ﬁrst exoplanet candidate astrometrically detected by Gaia to be successfully conﬁrmed based on RV follow-up observations.


Introduction
The publication of Gaia Data Release 3 (DR3) on June 13th, 2022 (Gaia Collaboration 2023a) has provided the first catalogue of >800 000 non-single star solutions in the astrometric, spectroscopic, and photometric channels (Gaia Collaboration 2023b).Of these, 169 227 are full, single-Keplerian astrometric orbital solutions (Gaia Collaboration 2023b; Halbwachs et al. 2023).The sample is dominated by solutions with inferred companion masses in the stellar regime, but about 1% of them (1915) corresponds to primaries orbited by objects in the sub-stellar mass regime (assuming zero flux ratio), with 72 companions with estimated masses formally <20 M Jup , that is, potentially planetary in nature (Gaia Collaboration 2023b).
The first sample of exoplanet candidates detected by Gaia astrometry contains a subset of nine previously known Dopplerdetected giant planets, which offered the means for their direct confirmation (Gaia Collaboration 2023b, and references therein).In a few cases, publicly available radial velocity (RV) data were used by Holl et al. (2023) to validate the orbital solu-Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).tions based on the close correspondence between the values of orbital period obtained from Gaia astrometry and the RV data.In this respect, the case of the companion orbiting HIP 66074 (Gaia DR3 1712614124767394816) constitutes a rather interesting conundrum.The Gaia solution has period P = 297 ± 2.8 d, eccentricity e = 0.46 ± 0.17, inclination i = 90 ± 5 • , and angular semi-major axis a 0 = 0.21 ± 0.03 mas.Using a reasonable guess to the mass of its K-dwarf primary, HIP 66074b is inferred to be a super-Jovian planet with M p ∼ 7.0 M Jup .Given the large companion mass, for an essentially edge-on orbit the expected RV semi-amplitude K from the Gaia solution would be K 300 m s −1 .In a work focussed on providing consistency tests between Gaia astrometry and Doppler data on selected systems with known and candidate exoplanets, Winn (2022) re-analysed the available Keck HIRES RVs of HIP 66074 (Butler et al. 2017).This work highlighted the presence of an obvious inconsistency: the RV-only solution has the same period of the Gaia orbit, but the K-value is 15 times smaller (∼20 m s −1 ).A combined analysis of the RVs and the Gaia solution for the source results in a very good fit, albeit requiring an unrealistically high flux ratio (Winn 2022).A more recent analysis performed by Marcussen & Albrecht (2023) achieved the same conclusions, ruling out the scenario of an astrophysical false positive produced by a binary system with a mass ratio almost equal to the flux ratio, based on the investigation of the public HIRES spectra.
Here, we present an analysis of precise RVs of HIP 66074 gathered at regular cadence with the HARPS-N spectrograph (Cosentino et al. 2012) at the Telescopio Nazionale Galileo (TNG) within the context of the programme Global Architecture of Planetary Systems (GAPS, Covino et al. 2013;Desidera et al. 2013).The combination of the new orbital solution and detailed simulations of Gaia observations of the system over the DR3 time-span, enabling characterisations of the existing biases in the Gaia-only solution, allows us to conclusively confirm the planetary nature of the candidate.Thus, we identified HIP 66074b/Gaia-3b1 as the first detected exoplanet by Gaia astrometry and the first-ever unambiguously confirmed astrometric detection of an exoplanet at any wavelength.

Spectroscopic observations
We planned an intensive RV monitoring of HIP 66074 with HARPS-N, collecting a total of 60 spectra between 4 June 2022 and 30 April 2023 (330 days), with a typical exposure time of 900 s.The spectra were reduced with version 3.7.1 of the HARPS-N Data Reduction Software (DRS) pipeline, which is maintained by the Italian centre for Astronomical Archive (IA2) 2 .We derived the RVs using v1.8 of the Template Enhanced Radial velocity Reanalysis Application (TERRA) pipeline (Anglada-Escudé & Butler 2012), particularly suited for late-type dwarfs such as HIP 66074 (Perger et al. 2017).The RV time series obtained with TERRA (reported in Table D.1) has a median of σ RV = 1.25 m s −1 and rms of 22.3 m s −1 .

Updated stellar parameters
Assuming as the input parameters the effective temperature (T eff ) and surface gravity (log g) from the StarHorse2 catalogue (Anders et al. 2022), we derived the final T eff , log g, and iron abundance ([Fe/H]) using both the spectral synthesis and equivalent width methods to the co-added spectrum of the target.In the first case, we considered the SME code (Piskunov & Valenti 2017;version 2020), while in the second case we considered the MOOG code (Sneden 1973;version 2019) and the iron line list by Biazzo et al. (2022).In both cases, we fixed the microturbulence velocity ξ to 0.5 km s −1 from the relationship by Adibekyan et al. (2012), and we used MARCS (Gustafsson et al. 2008) and ATLAS9-ODFNEW (Castelli & Kurucz 2003) grids of model atmospheres, obtaining consistent results.Moreover, within the spectral synthesis procedure, we adopted the macroturbulence velocity of 1.4 km s −1 by Brewer et al. (2016).We also derived T eff considering the line-depth ratio (LDR) method and appropriate calibrations LDR-T eff developed at the same resolution as HARPS-N (see Biazzo et al. 2011, and references therein).In all three cases, we obtained similar results within the uncertainties.Mean final values of the derived parameters are: T eff = 4300 ± 60 K, log g = 4.58 ± 0.06, and [Fe/H] = 0.12 ± 0.05 (see Table A.1).As a by-product, we also derived through the spectral synthesis technique the projected rotational velocity of v sin i = 1.8 ± 0.6 km s −1 .No lithium line was detected in the spectrum, thereby indicating the star is not young.
We derived the stellar mass, radius, and age based on a fit to the star's broad-band spectral energy distribution (SED) from the optical to the mid-infrared (see Fig. B.1) using the EXO-FASTv2 code (Eastman 2017;Eastman et al. 2019).The SED fit (see e.g.Stassun & Torres 2016) was performed using archival broad-band Tycho-2 and Johnson-Cousins B-and V-band magnitudes, i-band SDSS photometry, 2MASS JHK s near-infrared magnitudes, and WISE W1−W4 IR magnitudes.Within EXO-FASTv2 we utilised the YY-isochrones (Yi et al. 2001) to obtain the following stellar properties (also reported in Table A.1): M = 0.705 −0.023 −0.025 M , R = 0.690 +0.012 −0.013 R , and t = 7.9 +4.9 −4.1 Gyr.The space velocities support the membership to the thin disk and are well outside the kinematic space of young stars (Montes et al. 2001).These characteristics indicate a plausible age range of ∼2−8 Gyr, consistent with the X-ray non detection in the ROSAT All-Sky Survey.From the HARPS-N spectra we obtained a median Mount Wilson S -index value of S MW = 0.78 (the full S MW time-series is reported in Table D.1), corresponding (using the formalism described in Astudillo-Defru et al. 2017) to log(R HK ) = −4.80 ± 0.04 and an expected rotation period P rot = 31±3 d.This is close to the tentative rotation period of ∼35 d inferred from the analysis of the HARPS-N RV and S -index time series (Sect.3.2).Such a period is slightly longer than that of stars of similar colour in the 4 Gyr-old open cluster M 67 (Gruner et al. 2023), suggesting an older age.In summary, all the indicators consistently support an old age.

Spectroscopic orbital solution
We initially fitted a Keplerian orbit to the HARPS-N RV measurements using the publicly available Monte Carlo (MC) nested sampler and Bayesian inference tool MULTINEST V3.10 (e.g.Feroz et al. 2019), through the pyMULTINEST wrapper (Buchner et al. 2014) We performed the same analysis on the S -index time-series, which also returned a clear periodicity at ∼35 d (Fig. D.1, bottom panel).We also downloaded the publicly available photometric data of HIP 66074 from the Mikulski Archive for Space Telescopes (MAST) portal3 that were gathered by the TESS mission in nine non-consecutive sectors and performed a periodogram analysis of the light curves extracted with an independent pipeline (see e.g.Nardiello et al. 2022).Although TESS observations are not designed to easily detect P rot > 15 d, we found evidence of rotational modulation in the range 32−41 d upon inspection of the TESS photometry of different sector subsets, or taken as a whole (Fig.

B.2).
Given the convergent view on the likely presence of rotational modulation in the HARPS-N RVs, we then expanded the model to include a stellar activity term using the Gaussian process (GP) regression package GEORGE (Ambikasaran et al. 2015) and adopting a quasi-periodic (QP) kernel (see, e.g.−0.04 Damasso et al. 2020).We included the publicly available HIRES RVs (omitting, as in Winn 2022, the spectrum gathered at epoch JD = 2455042.76395) in the modelling of the Keplerian signal (requiring the addition of an RV offset γ HIRES and an uncorrelated jitter σ jit HIRES to the set of model parameters); however, the GP QP model was applied only to the HARPS-N RVs, given the sparseness of the former dataset.The combined HARPS-N+HIRES time baseline exceeds 14 years and it is therefore suitable to look for evidence of any long-term RV trends possibly due to long-period companions consistently present in both datasets.In addition to the single Keplerian + GP QP model, we then tested for the presence of acceleration and curvature in the data, but to no avail.The solution with single Keplerian + GP QP is favoured in terms of Bayesian evidence, with marginal likelihood differences of ∆ ln Z > +4, and we adopted it as fiducial.The priors and final results of the global modelling of the HARPS-N RVs are reported in  This result clearly helps resolving the puzzling discrepancy with the Gaia orbital solution initially highlighted in the literature, but not entirely.Based on the solid evidence obtained with the HARPS-N RVs, we next investigate the specific biases present in the Gaia-only orbital solution.

Understanding biases in the Gaia-only solution
In order to fully reconcile the still-existing discrepancy between our RV-measured orbit and the published Gaia DR3 solution, we performed a set of numerical simulations.We generated a set of 140 synthetic Gaia astrometric observations of HIP 66074.The transit times (determined relative to the Gaia DR3 reference epoch T ref = 2016.0),scan angles, and along-scan parallax factors of the Gaia observations of HIP 66074 over the DR3 time-span were obtained from the Gaia Observation Forecast Tool (GOST) 4 .Each of the synthetic time-series of along-scan coordinates w was produced with the stellar motion described by the five astrometric parameters from Table A.1.The astrometric signal induced by a single companion was then added linearly, using the stellar mass value as determined in Sect.3.1 and the orbital parameters P, T 0 , e, and ω and the derived minimum mass from Sect.3.2, varied within their respective uncertainties.The two remaining orbital elements, the orbital inclination i and the longitude of the ascending node Ω, were generated from uniform distributions over the [0, π] range, and we then computed the true companion mass and angular semi-major axis based on i.Finally, the w measurements were perturbed by Gaussian random uncertainties with a standard deviation σ w = 0.1 mas, appropriate for a star of magnitude similar to HIP 66074 for DR3-level astrometry (see Holl et al. 2023).Each of the 140 astrometric time series was then fitted with a 12-parameter model (five astrometric parameters, seven orbital elements) using a Markov chain Monte Carlo (MCMC) algorithm, adopting the emcee Affine Invariant MCMC Ensemble sampler by Foreman-Mackey et al. (2013).We adopted a partly linearised Keplerian orbit model (e.g.Holl et al. 2023;Halbwachs et al. 2023).For the non-linear parameters T 0 , P, e, we used uniform priors to avoid biasing the orbital solution, as illustrated in Table C.1.
The results are illustrated in the four panels of Fig. 2. As we can see, the medians and standard deviations of the fitted eccentricity and derived inclination are e fit = 0.53 ± 0.13 and i fit = 93 ± 27 • .The combination of the intrinsically onedimensional Gaia w measurements and small size of the astrometric perturbation induces a systematic underestimation of the eccentricity (an issue already hinted at by Holl et al. 2023) and a bias against close-to-face-on orbits.Furthermore, the retrieved semi-major axis is systematically overestimated, with fitted semi-major axis a fit σ w , except for very small i values corresponding to simulated a σ w .For quasi-face-on configurations (e.g.simulated i 5 • and i 175 • for a retrograde and prograde orbits, respectively), the fitted eccentricity gets closer to the input value, but the pronounced suppression of quasi-faceon solutions remains.This bias is due to the partially linearised model fitted to noisy data, which it has already been qualitatively discussed by Gaia Collaboration (2023b).Finally, we notice that the median significance of a fit /σ a when a σ w is typically around 3, well below the value of ∼7 of the published Gaia solution.In our simulations, this level of significance is obtained for a σ w , which implies that the size of the published Gaia orbit is possibly correct, but it could also be overestimated by up to a factor of ∼2.
We then go on to consider that the semi-major axis of the primary (in mas) obtained by Gaia can be written as (Pourbaix 2001): where K is in m s −1 , P in yr, and in mas.As the orbital elements from the spectroscopic orbit (including the K-value) are now robustly determined, we can make a better-informed statement on the true inclination of the orbit.If a 0 = 0.21 ± 0.03 mas derived from the Gaia solution is correct, then this implies i 6.5 • .If a 0 is overestimated by as much as a factor of 2, then i 13.0 • .We cannot completely rule out the possibility of a significantly higher inclination and, therefore, a true mass closer to the minimum mass.However, this would imply (see Fig. 2) a in the range 60−25 µas for i between 25 • and 90 • , respectively, significantly below the quoted measurement uncertainties.It is unclear whether such a small perturbation size would have been effectively detectable in Gaia DR3 astrometry.We therefore conclude that the true orbit orientation indeed corresponds to a close to face-on configuration and the true mass estimate of HIP 66074b/Gaia-3b is likely in the range 3−7 M Jup (close to the one derived by Gaia DR3 astrometry), fully resolving the tension between the Gaia solution and Doppler spectroscopy.

System architecture
HIP 66074 is known to have a faint companion, 2MASS J13324530+7459441, at an angular separation of 43.8 , as first discussed by Gomes et al. (2013).The Gaia DR3 astrometry firmly confirms the physical association of the two objects, albeit with a statistically significant (∼15σ) proper motion difference L15, page 4 of 10 ∆µ in right ascension.Gomes et al. (2013) derived an L2 spectral type from optical spectroscopy, T eff = 2080 ± 260 K, and bolometric luminosity log L/L = −3.78± 0.045.We took advantage of the stellar characterisation of the primary (Sect.3.1) for further inferences.From the adopted stellar age and the models by Baraffe et al. (2015), we obtained a mass of 0.076±0.001M for the companion (considering the Gaia and near-infrared absolute magnitudes and the spectroscopic T eff ; the uncertainty does not include the systematic errors of the models and the small effect of non-solar metallicity), just above the Hydrogen-burning limit.We then concluded that HIP 66074 has a very low-mass stellar companion at a projected separation of ∼1550 au.
The existence of a wide companion around a star with a planet with extreme eccentricity cannot be considered a surprise.Indeed, all the previously known planetary companions with eccentricities larger than 0.9 have further companions on wide orbits (HD 4113, Tamuz et al. 2008;HD 7449, Cheetham et al. 2018;HD 80606, Naef et al. 2001;HD 20782, Desidera & Barbieri 2007), and, more generally, a link between binary fraction and planet eccentricity has been reported (e.g.Moutou et al. 2017;Su et al. 2021).This points to a key role of dynamical interactions with an outer perturber in the generation of extreme planet eccentricities (Mustill et al. 2022).To further detail the mechanism of such interactions, we derived the timescale for the Kozai modulation, following Takeda & Rasio (2005).This is shown to be comparable to the age of the system for very high eccentricities of the outer orbit ( 5Gyr for e = 0.9) and longer than the age of the universe at lower eccentricities (25 Gyr for e = 0.67).We computed γ, the angle between the separation vector and the relative velocity vector of the binary (Tokovinin & Kiyaeva 2016), finding γ = 67 ± 4 • .This indicates, in a statistical sense, a not particularly high binary eccentricity (e.g.Tokovinin & Kiyaeva 2016;Hwang et al. 2022).Furthermore, we verified both analytically (Misner et al. 1973) and via direct numerical integration that the relativistic precession timescale of the orbit is much shorter than that of the Kozai modulation, which therefore is effectively inhibited.The wideseparation low-mass stellar companion of HIP 66074 is thus unlikely to be the one responsible for the observed eccentricity of Gaia-3b.Other scenarios such as planet-planet scattering events (Carrera et al. 2019) might be more realistic scenarios for explaining this discrepancy.

Conclusions
We confirm the first Gaia astrometric planet discovery, HIP 66074b/Gaia-3b, based on new HARPS-N RV data.The object was at the centre of a much-debated, puzzling discrepancy between Gaia DR3 astrometry and Keck HIRES RVs of HIP 66074.Both the Gaia solution and HIRES RVs indicate the presence of a planetary-mass companion with a period of ∼300 days, but with entirely incompatible measurements of the astrometric orbit size and RV semi-amplitude, given the remainder of the orbital elements.This conundrum was fully resolved based solely on the critical contribution of a dense RV-monitoring with HARPS-N RV that allowed for a sampling of the periastron passage of the second-highest eccentricity gas giant ever found.In turn, based on detailed numerical simulations, this allowed the discrepancy with the Gaia-only solution to be reconciled.The Gaia orbit suffers from two main biases: the eccentricity is underestimated and the edge-on configuration is incorrect, while the true inclination angle is small, namely, i 13 • .The angular orbit size is likely correct or overestimated by a factor of at most ∼2.The true mass estimate of Gaia-3b, between approximately 3 and 7 M Jup , unambiguously identifies the object as a super-Jovian planet, the first-ever astrometrically detected companion fully in the planetary regime 5 to be independently confirmed by another technique after many decades of attempts (e.g.Sozzetti & de Bruijne 2018, and references therein).
Gaia-3b joins the small sample of giant exoplanets with extremely high eccentricities (e > 0.9), with the peculiarity that it has been found to be orbiting the star with the lowest mass star of the lot.The clearly identified wide, very-low-mass stellar companion, straddling the threshold between Deuterium-and Hydrogen-burning objects, does not appear to effectively influence the orbit of Gaia-3b via Kozai cycles, while Jovian-mass (or larger) companions out to 5−10 au are also ruled out by the lack of detectable long-term RV trends and statistically significant Hipparcos-Gaia DR3 proper motion anomaly (Kervella et al. 2022;Brandt 2021).The HIP 66074 planetary system therefore constitutes an excellent laboratory for in-depth studies on the influence of a third body and of the effects of tidal circularisation on the parameters of the exoplanet Gaia-3b.Further investigations of the system are warranted and these will benefit from the effective combination of RVs and Gaia time-series astrometry spanning 5.5 years, when the next major data release, DR4, is published near the end of 2025.
, with uniform priors on the model parameters.The orbital model has the following free parameters: the epoch of inferior conjunction T 0, b ; the orbital period P b ; a residual RV offset γ HN ; the RV semi-amplitude K b ; √ e b cos ω , b and √ e b sin ω , b , where e b is the eccentricity and ω , b the argument of periastron; an uncorrelated jitter term σ jit, HN added in quadrature to the formal uncertainties.We obtained K ∼ 94 m s −1 , P ∼ 285 d, and e ∼ 0.94.A Generalized Lomb-Scargle (GLS; Zechmeister & Kürster 2009) periodogram analysis of the residuals showed the presence of a peak at ∼35 d, albeit with a high bootstrap-based false alarm probability of ∼11% (Fig. D.1).

Fig. 1 .
Fig. 1.Phase-folded HIRES and HARPS-N RVs superposed to the bestfit Keplerian orbit (blue curve), calculated using the median values of the posteriors.Phase zero (or phase one) corresponds to the time of inferior conjunction.The residuals of the best-fit model are shown in the bottom panel.

Fig. 2 .
Fig. 2. Gaia astrometry simulations results.From top to bottom: (a) derived inclination vs. simulated value; (b) fitted eccentricity vs. simulated value; (c) derived (blue) and simulated (red) astrometric semimajor axis vs. injected orbital inclination; the black dotted line shows the semi-major axis of the Gaia solution, a 0 ; (d) significance of the a fit vs. simulated value; the black dotted line indicates the value obtained for the Gaia solution.In all panels yellow diamonds highlight orbits with inclination within 5 • of an exactly face-on configuration.
TableD.1 lists the TERRA-extracted HARPS-N RV data and their uncertainties, along with the time series of the S MW spectroscopic activity indicator.Figure D.1 shows GLS periodograms of the HARPS-N RV residuals to a single Keplerian fit and of the S -index time-series.

Fig
Fig. D.1.GLS periodogram of the residuals to a Keplerian fit of the HARPS-N RVs alone (top).GLS periodogram of the S MW time-series (bottom).

Fig. E. 2 .
Fig. E.2.Posterior distributions of the model (hyper)parameters of our assumed best-fit model, including a Keplerian for planet b and a GP quasi-periodic correlated stellar activity signal fitted to the HARPS-N RVs of HIP 66074.

Table 1 .
Global modelling of the HIP 66074 HIRES and HARPS-N RVs: priors, best-fit, and derived parameters.
Table 1, whileFig. 1 shows the best-fit Keplerian solution overplotted to the phase-folded RVs and Figs.E.1 and E.2 show the best-fit solution for the stellar activity component and the full set of joint posteriors, respectively.An inspection of Fig. E.2 shows multiple peaks in the O'Toole et al. 2009ion for P. We should keep in mind that the very sparse distribution of HIRES RVs of HIP 66074 (on average, two data points distributed over a full orbital period, never sampling phases close to periastron) is prone to the introduction of period aliases particularly in the case of such a high-eccentricity orbit (see e.g.O'Toole et al. 2009).The results reported in Table1seem to indicate coherence of the activity signal (large corre-

Table D .
1. HARPS-N RVs and S MW time-series of HIP 66074.