# Astrometric detection of a Neptune-mass candidate planet in the nearest M-dwarf binary system GJ65 with VLTI/GRAVITY

R. Abuter, A. Amorim, M. Benisty, J. P. Berger, H. Bonnet, G. Bourdarot, P.

Bourget, W. Brandner, Y. Clénet, R. Davies, et al.

► To cite this version:

R. Abuter, A. Amorim, M. Benisty, J. P. Berger, H. Bonnet, et al. Astrometric detection of a Neptune-mass candidate planet in the nearest M-dwarf binary system GJ65 with VLTI/GRAVITY. Astronomy & Astrophysics - A&A, 2024, 685, 10.1051/0004-6361/20249547. insu-04836871

HAL Id: insu-04836871

https://insu.hal.science/insu-04836871v1

Submitted on 13 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
LETTER TO THE EDITOR

# Astrometric detection of a Neptune-mass candidate planet in the nearest M-dwarf binary system GJ65 with VLTI/GRAVITY

GRAVITY Collaboration*: R. Abuter⁴, A. Amorim⁸,⁷, M. Benisty³, J. P. Berger³, H. Bonnet⁴, G. Bourdarot¹,²,*, P. Bourget⁴, W. Brandner⁵, Y. Clénet², R. Davies¹, F. Delplancke-Ströbele⁴, R. Dembet², A. Drescher¹, A. Eckart⁶,¹⁴, F. Eisenhauer¹,¹⁶, H. Feuchtgruber¹, G. Finger¹, N. M. Förster Schreiber¹, P. Garcia¹⁰,⁷, R. Garcia-Lopez¹⁷,¹⁸, F. Gao¹¹, E. Gendron², R. Genzel¹,¹², S. Gillessen¹, M. Hartl¹, X. Haubois⁹, F. Haussmann¹, T. Henning⁵, S. Hippler⁵, M. Horrobin⁶, L. Jochum⁹, L. Jocou³, A. Kaufer⁹, P. Kervella²,²*, S. Lacour², V. Lapeyrère², J.-B. Le Bouquin³, C. Ledoux⁹, P. Léna², D. Lutz¹, F. Mang¹, A. Mérand⁴, N. More¹, M. Nowak¹⁹, T. Ott¹, T. Paumard², K. Perraut³, G. Perrin², O. Pfuhl⁴,²*, S. Rabien¹, D. C. Ribeiro¹, M. Sadun Bordoni¹, J. Shangguan¹, T. Shimizu¹, J. Stadler¹³, O. Straub¹⁵, C. Straubmeier⁶, E. Sturm¹, L. J. Tacconi¹, K. R. W. Tristram⁹, F. Vincent², S. von Fellenberg¹⁴, F. Widmann¹, E. Wieprecht¹, J. Woillez⁴, S. Yazici¹, and G. Zins⁴

*(Affiliations can be found after the references)*

Received 8 February 2024 / Accepted 20 March 2024

## ABSTRACT

The detection of low-mass planets orbiting the nearest stars is a central stake of exoplanetary science, as they can be directly characterized much more easily than their distant counterparts. Here, we present the results of our long-term astrometric observations of the nearest binary M-dwarf Gliese 65 AB (GJ65), located at a distance of only 2.67 pc. We monitored the relative astrometry of the two components from 2016 to 2023 with the VLTI/GRAVITY interferometric instrument. We derived highly accurate orbital parameters for the stellar system, along with the dynamical masses of the two red dwarfs. The GRAVITY measurements exhibit a mean accuracy per epoch of 50–60 ms in 1.5 h of observing time using the 1.8 m Auxiliary Telescopes. The residuals of the two-body orbital fit enable us to search for the presence of companions orbiting one of the two stars (S-type orbit) through the reflex motion they imprint on the differential A–B astrometry. We detected a Neptune-mass candidate companion with an orbital period of *p* = 156 ± 1 d and a mass of *m*<sub>p</sub> = 36 ± 7 M<sub>⊕</sub>. The best-fit orbit is within the dynamical stability region of the stellar pair. It has a low eccentricity, *e* = 0.1–0.3, and the planetary orbit plane has a moderate-to-high inclination of *i* > 30° with respect to the stellar pair, with further observations required to confirm these values. These observations demonstrate the capability of interferometric astrometry to reach microarcsecond accuracy in the narrow-angle regime for planet detection by reflex motion from the ground. This capability offers new perspectives and potential synergies with Gaia in the pursuit of low-mass exoplanets in the solar neighborhood.

**Key words.** instrumentation: high angular resolution – instrumentation: interferometers – astrometry – planets and satellites: detection – planets and satellites: dynamical evolution and stability – stars: low-mass

## 1. Introduction

M-dwarfs are favorable targets in the search for signatures of planets thanks to their low stellar masses (0.1–0.25 M<sub>☉</sub>), as the presence of planets induces a significantly larger reflex motion on the parent star than for sun-like stars. This makes them ideal targets both for radial velocity and astrometric searches. In our close solar neighborhood, M-dwarfs represent 80% of the stellar population (Winters et al. 2019; Chabrier & Baraffe 2000). Thus, the characterization of the properties of the stellar host and their monitoring are of key interest. In that respect, radial

velocity and astrometry probe complementary parameter spaces. Radial velocity (RV) measurements are most sensitive to close-in planets, but are ultimately limited to a precision of a few m s⁻¹ in M-dwarfs due to stellar variability and flaring events (Meunier 2021). Conversely, astrometry probes planetary companions at larger separation (typically 0.1 to a few astronomical units, au) and is less sensitive to flare activity, with a typical noise contribution of 10 µau for M-type stars (Eriksson & Lindegren 2007; Lagrange et al. 2011). In addition, in the case of M stars, the flux contribution of flares decreases by two orders of magnitude in the K band (2.0 µm–2.4 µm) compared to the UV and optical (Davenport et al. 2012), favoring the infrared domain for precision astrometry for these objects. The typical astrometric signature of a planet with a few Earth-mass on a 0.1 au orbit around an M dwarf located at 1 pc is on the order of 1–10 µas, below the present limit of space and ground based astrometric instrumentation.

In this context, narrow-angle astrometry with interferometry (Eisenhauer et al. 2023) provides a powerful technique for

<footnote>* GRAVITY is developed in a collaboration by MPE, LESIA of Paris Observatory/CNRS/Sorbonne Université/Univ. Paris Diderot and IPAG of Université Grenoble Alpes/CNRS, MPIA, Univ. of Cologne, CEN-TRA – Centro de Astrofísica e Gravitação, and ESO.</footnote>

<footnote>** Corresponding authors: G. Bourdarot (e-mail: bourdarot@mpe.mpg.de), P. Kervella (e-mail: pierre.kervella@obspm.fr), and O. Pfuhl (e-mail: opfuhl@eso.org).</footnote>
monitoring stellar orbits with very high precision. Shao & Colavita (1992) first discussed the potential of micro-arcsecond precision from the ground with long-baseline interferometry to observe two objects contained in the same isopistonic area. This technique was first demonstrated at PTI (Colavita et al. 1999), achieving 100 µas accuracy (Muterspaugh et al. 2010), but put on hold after the cancellation of the large space-based NASA mission SIM (Unwin et al. 2008) and first-generation ground-based instruments (Keck/ASTRA, Wouillez et al. 2010; VLT/PRIMA, Delplancke et al. 2006). The full potential of dual-field interferometry was finally reached with GRAVITY (GRAVITY Collaboration 2017), routinely achieving an accuracy of a few 10 µas on 8 m class telescopes for objects as faint as $m_K = 19.5$ mag (GRAVITY Collaboration 2018, 2022a).

In the following, we present the astrometric monitoring of GJ65 AB (also known as: Luyten 726–8, UV Ceti+BL Ceti, WDS J01388–1758 AB) we conducted with VLTI/GRAVITY on the 1.8 m VLTI Auxiliary Telescopes between 2016 and 2023. GJ65 AB is a binary star located at $d = 2.67$ pc, and one of the Sun’s nearest neighbors. GJ65 is the class prototype of UV Ceti Type stars (Luyten 1949; Joy & Humason 1949). The separation of the binary ($a = 5.4$ au) is wide enough to neglect magnetic and tidal interactions between the stars, but its orbital period $P = 26.38$ yr is short enough to allow for an extremely accurate determination of the orbital parameters (Kervella et al. 2016). This makes GJ65 a prime target in determining the masses of the two stellar components. The stars’ properties are very similar to Proxima Centauri (Kervella et al. 2016). The components A and B in GJ65 have nearly identical masses (≈0.12 $M_\odot$). No RV detection of a planet has been reported, due to the fast rotation of the two stars and their strong variability through flaring events that limit the accuracy of RV measurements. In astrometry, indications of the presence of a third body from the proper motion anomaly were found from Gaia DR2 and EDR3 but with no clear planet signature (Kervella et al. 2019, 2022). Understanding planet formation in close binary systems is still challenging and models show that it is hindered by the dynamical excitation and gravitational instability (Holman & Wiegert 1999; Marzari & Thebault 2019). Therefore, obtaining accurate measurements of the stellar orbit and constraining the presence of potential companions is imperative to testing these models.

The outline of the Letter is as follows. Section 2 describes the observations of our monitoring program and Sect. 3 presents the principle of dual-field observations used for narrow-angle astrometry with long baseline interferometry, and the performance achieved in the program. In Sect. 4, we present the orbital parameters and the analysis of the residuals of the astrometry of GJ65. Finally in Sect. 5, we discuss the properties and the dynamical stability of this system, and offer an outlook for future studies.

## 2. Observations

We observed GJ65 from 2016 to 2023 with GRAVITY/ATs, with typically 15 h of observations per semester. The individual epochs have a total observing time ranging from 1 h to 3.3 h including overheads. All data were acquired in medium resolution ($R = 500$) and SPLIT polarization mode to separate the stellar light in two linear polarizations. Each observation consists of a series of field swaps of about 20 min each, during which the two fibers observe each star successively in order to extract the differential astrometric signal. The principle of measurement in dual-field observations and of swap measurement is detailed in Sect. 3. The observation sequence does not require an interferometric transfer function calibrator given that the visibility amplitudes are not used in dual-field astrometry and the phases are self-calibrated by the swapping procedure.

## 3. Phase-referenced astrometry

Narrow-angle astrometry is enabled by the dual-field design of GRAVITY (GRAVITY Collaboration 2017). In this mode, the interferometric beam combiner is able to simultaneously observe two separated interferometric field of view, with two fibers for each telescope, thereby providing two interference patterns for each baseline. The angular separation of the two objects is measured from the optical path difference (OPD) between the white light fringe of these two fringe patterns, delivering astrometry with a high level of precision that is defined by the interferometric baseline. In order to compensate for additional perturbations on the path between these two channels, a metrology system was implemented (Gillessen et al. 2012; Lippa et al. 2016) that measures the OPD from the two-beam combiner to the top of the telescopes. The relevant baseline in astrometric observations is the narrow angle baseline (Wouillez & Lacour 2013; Lacour et al. 2014), defined by the end of the metrology point, located on the secondary mirror spider arms of each telescope. The relation between the phase measured by GRAVITY and the angular separation between the targets is then (Lacour et al. 2014):

$$ \phi = (\phi_{\text{SC}} - \phi_{\text{FT}}) - \phi_{\text{MET}} - \frac{2\pi}{\lambda} D(\lambda_{\text{m}}, \lambda) \\ = \frac{2\pi}{\lambda} (s_{\text{B}} - s_{\text{A}}) \cdot \mathbf{B}_{\text{NAB}} + \phi_{0}(\lambda) \quad (1) $$

with $\phi_{\text{SC}}$ as the phase measured on the science channel of GRAVITY, $\phi_{\text{FT}}$ the phase measured on the fringe-tracker, $\phi_{\text{MET}}$ the differential OPD measured by the metrology, $\mathbf{B}_{\text{NAB}}$ the narrow-angle baseline, $D(\lambda_{\text{m}}, \lambda_{\text{SC}})$ the dispersion between the wavelength of the metrology and of the starlight, $s_{\text{A}}$ and $s_{\text{B}}$ the position vector of stars A and B in the plane of sky, $\lambda$ the wavelength of the science channel and $\phi_{0}(\lambda)$ the zero point of the metrology.

To carry out the phase astrometry, it is necessary to remove $\phi_0(\lambda)$, the constant internal offset of the metrology. In dual-field observations, this can be done either by choosing a reference point to the observations, which will be used as phase reference and subtracted from the astrometry (GRAVITY Collaboration 2018, 2019) or by measuring this zero-point independently by the means of successive swaps between the two fibers. The latter option is used in the case of GJ65. We detail the methodology of the measurement in Appendix A.

We implemented this method and applied it on the standard reduced files produced by the GRAVITY reduction pipeline (Lapeyrere et al. 2014). The output of the process is finally the relative angular separation $s_{\text{B}} - s_{\text{A}}$ of the A and B components.

## 4. Analysis

### 4.1. Stellar orbit

We started by fitting the orbital elements of the stellar binary based on the astrometric data. The orbital fit includes our new GRAVITY data, complemented by archival astrometry data from WDS catalog and NACO published in Kervella et al. (2016). The result of the fit is dominated by the GRAVITY data, as shown in Fig. 1 given the much higher accuracy on these points, roughly a factor 100 compared to the archival data. For each GRAVITY point, the uncertainty and the degree of correlation of the error
<figure><img src="image_2.png" /><figcaption>Fig. 1. Left: Orbit of GJ65 AB, combining WDS, NACO and GRAVITY data. Right: residuals of the stellar orbit.</figcaption></figure>

in RA and Dec are provided as part of the GRAVITY data reduction. The covariance matrix can be reconstructed from these quantities and is integrated in the orbital fit. The result of the fit is shown in Table 1. We expressed the orbital parameters following the same convention as used in Gaia Collaboration (2023), Halbwachs et al. (2023), and shown in Appendix B. The dynamical mass of the binary is constrained by the combination of the GRAVITY astrometry with UVES radial velocity retrieved from ESO public archive, whose joint fit with the astrometry is shown on Figs. D.2 and D.3. This study refines the initial estimate without dual-field astrometry data given in Kervella et al. (2016), with which they agree. We note the discrepancy in the value of the inclination in Kervella et al. (2016) due to the difference of convention used in our current Letter, which is the same as in Gaia DR3.

## 4.2. Periodogram

We computed a generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009) on the residuals of the astrometry of the stellar orbit. The periodograms of the residuals in right ascension (RA) and declination (Dec) were first computed independently, taking into account their error bars. In addition, we computed the joint periodogram of the signal following the methodology described by Catanzarite et al. (2006), which consists of summing the periodograms in RA and Dec, however the detection is mainly driven by the signal in Dec. The result is shown in Fig. 2. We used a bootstrapping approach to compute the false alarm probability (FAP) and the detection threshold. To do so, a large number ($N = 10^4$) of random periodograms in RA and Dec was generated, then drawn with replacement at the measurement epochs of our sample and following a normal distribution with a zero mean and a standard deviation corresponding to the $1\sigma$ uncertainty of that point. For each draw, a periodogram

<table><caption>Table 1. Orbital parameters of GJ65 A and B obtained with GRAVITY.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">Orbital parameters:</td></tr><tr><td>Period <i>P</i> (yr)</td><td>26.38 ± 0.002</td></tr><tr><td><i>a</i> (AU)</td><td>5.459 ± 0.002</td></tr><tr><td><i>i</i> (deg)</td><td>128.0 ± 0.1</td></tr><tr><td><i>e</i></td><td>0.6172 ± 0.0001</td></tr><tr><td>Ω (deg)</td><td>325.9 ± 0.1</td></tr><tr><td>ω (deg)</td><td>103.2 ± 0.1</td></tr><tr><td><i>T</i><sub>0</sub> (MJD)</td><td>41333 ± 8</td></tr><tr><td colspan="2">Masses:</td></tr><tr><td><i>m</i><sub>A</sub> (<i>M</i><sub>☉</sub>)</td><td>0.122 ± 0.002</td></tr><tr><td><i>m</i><sub>B</sub> (<i>M</i><sub>☉</sub>)</td><td>0.116 ± 0.002</td></tr><tr><td><i>m</i><sub>A</sub> + <i>m</i><sub>B</sub> (<i>M</i><sub>☉</sub>)</td><td>0.238 ± 0.003</td></tr><tr><td>Parallax <sup>(1)</sup> (mas)</td><td>371.92 ± 0.42</td></tr></tbody></table>

Notes. <sup>(1)</sup>From Kervella et al. (2022).

was generated in RA and Dec. The FAPs were finally estimated by computing the 10%, 1%, and 0.15% quantiles of the set of generated periodograms.

Finally, in order to distinguish the impact of the temporal sampling on the final periodogram, we computed the periodogram of the sample replaced by constant values, superimposed on Fig. 2 (gray curve). This “blank” periodogram shows the typical frequencies of our sampling, with a distinct peak at $p = 365$ d for the one-year annual periodicity of the program.

The periodogram shows a peak at $p = 156$ d above a FAP = 0.15% threshold. The peak does not overlap with the test periodogram generated on white noise at the temporal sampling of our data set, nor at a harmonic of the one-year period associated to the main period of our program (see Fig. 2). For the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Generalized Lomb-Scargle periodogram. Top: periodogram in RA (orange) and Dec (blue) computed for the residuals of the binary orbit. Bottom: residuals after the subtraction of the planet candidate at *p* = 156 d.</figcaption></figure>

moment, we have discarded the peaks with a period <30 d, for which we consider our temporal sampling is not sufficient.

### 4.3. Orbit fitting

In a second step, we performed an orbital fit of the residuals to obtain a second detection of candidate, using a methodology that is different from the Lomb-Scargle periodogram, and to determine the orbital parameters of the companion planet. Our astrometric model describes the reflex motion of a companion in Keplerian motion around one of the components of the system. With respect to both our stellar orbit and the planetary orbit, the four orbital elements (*a*, *i*, Ω, ω) can be replaced by a linear model using the Thiele-Innes elements.

Firstly, we performed an orbital fit over a grid of fixed periods and fit all the other orbital parameters. In this way, we fixed the period, which is the most non-linear parameters of the fit. Since *P* is fixed, the fit is left with only two non-linear parameter (*e*, *t*₀) to solve the Kepler equation, and the four Thiele-Innes parameters (A, F, B, G), which are linear coefficients of the elliptical rectangular coordinates (*X*, *Y*). This method allows us to make the fit robust by limiting the number of non-linear parameters. The covariance matrix of the GRAVITY data, assuming the noise is random Gaussian, is taken into account in the fit as:

$$ \chi^2 = (\mathbf{y} - \mathbf{y}_{\text{model}})^T \Sigma^{-1} (\mathbf{y} - \mathbf{y}_{\text{model}}). \quad (2) $$

Finally, the Campbell parameters were deduced from the Thiele-Innes parameters, for which we again used the same convention as in Gaia Collaboration (2023), Halbwachs et al. (2023). The result of the grid search over periods converges to the same period *p* = 156 d as the Lomb-Scargle periodogram, as shown in Fig. B.1. The final orbital solution is obtained by allowing all parameters to vary, and converges to a Neptune-mass

planet with mass 35–40 $M_⊕$. The result for the best period and the phase folded orbit are shown on Fig. 3. The orbital parameters of the planet are shown in Table 2.

We note that the fit of the orbit is degenerate to four solutions in total: two possible host stars, and two possible longitudes of the ascending nodes (Ω, Ω + 180°) of the planet. However, the prograde solutions (two out of the four in total) are favored from a dynamical point of view, and are the ones shown in Table 2. The degeneracy on the host star could be lifted by measuring the absolute astrometry of the binary, and the degeneracy of the longitude of the ascending nodes by measuring the radial velocity of the host star caused by the motion of the planet. However, we note that the location of the planet in the plane of sky is known within two possible solutions (degeneracy of the host star), the degeneracy of the ascending node does not affect the knowledge of the location of the planet in the plane of sky (only *z*). The projection of the GJ65 system onto the sky is shown in Fig. D.1.

### 4.4. Additional companions and complementary RV data

Our orbital solution would lead to a typical RV amplitude of ≈20 m s⁻<fnref n="1" />. We used UVES data to constrain the stellar orbit (see Sect. 4.1) with a typical accuracy of the order of ±100 m s⁻<fnref n="1" />; however this is not sufficient to constrain a planetary companion. The two components of GJ65 are difficult targets for high precision RV measurements. Firstly, their very strong magnetic activity and fast rotation results in spotted and time-variable apparent surfaces (Barnes et al. 2017) that affect the spectral line profiles. Secondly, their relatively small angular separation (Fig. 1) induces crosstalk between the two stars for ground-based, seeing-limited spectrographs.

We note the presence of peaks at shorter periods ≈20 d both in the Lomb-Scargle periodograms and the period search. Such an astrometric signal would correspond to a Saturn-mass planet on a <0.1 AU orbit. Here, we limited our analysis to the periods >30 d due to the temporal sampling of our monitoring under one month. However, the typical RV amplitude of this signal is ≈100 m s⁻<fnref n="1" />, which could be confirmed by RV data. As a follow-up to the <30 d and the 156 d candidates, additional RV observations would provide complementary information in this system. The NIRPS high-stability spectrograph at La Silla (Bouchy et al. 2017) could be a prime instrument for GJ65 AB, thanks to its adaptive optics system which allows for the flux to be isolated from each components of the binary, and in the longer term the high-precision ELT/ANDES instrument.

## 5. Discussion

The candidate exoplanet that we detect in GJ65 is the second nearest after Proxima Centauri's exoplanet system at the date of writing<fnref n="1" />. It thus offers unique opportunities for the further characterization of this exoplanet. The astrometric detection allows for the mass, inclination and orbital parameters of the candidate to be derived. The location of the planet can be compared with the expected orbital stability region of test particles as derived by Holman & Wiegert (1999, see Eq. (1) of that paper; Appendix C). The orbital parameters of GJ65AB in Table 1 allow us to derive a critical semi-major axis of $a_c = 0.46$ au. The exoplanet candidate's orbit lies within this region, with a semi-major axis $a \approx 0.283$ au for the largest value (Table 2). This candidate could thus have formed and evolved within the region of

<footnote n="1">¹ exoplanetarchive.ipac.caltech.edu</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Orbital solution for the best period, *p* = 156 d (left). Phase-folded orbit (right) with the average data displayed as coloured dots.</figcaption></figure>

<table><caption>**Table 2.** Orbital parameters of the planet for the best period (reduced chi-square $\chi^2_r = 1.03$).</caption><thead><tr><th>Parameter</th><th>Star B</th><th>Star A</th></tr></thead><tbody><tr><td>Period (day)</td><td>156 ± 1</td><td>156 ± 1</td></tr><tr><td>Mass (M<sub>⊕</sub>)</td><td>36 ± 6</td><td>39 ± 7</td></tr><tr><td>a (AU)</td><td>0.274 ± 0.002</td><td>0.283 ± 0.002</td></tr><tr><td>e</td><td>0.27 ± 0.21</td><td>0.33 ± 0.30</td></tr><tr><td>i (deg)</td><td>89 ± 9</td><td>88 ± 6</td></tr><tr><td>Ω (deg)</td><td>9(*) ± 7</td><td>8(*) ± 7</td></tr><tr><td>ω (deg)</td><td>230 ± 36</td><td>248 ± 30</td></tr></tbody></table>

**Notes.** The host star is indicated for both solutions. (*)Both values (Ω, Ω + 180°) are compatible with this parameter given that it is degenerated (see text).

stability. The formation of planets at relatively short separation binaries (<10 au) is challenging due the dynamical interactions in these systems and the tidal truncation of the protoplanetary disks in early phases of formation, as a disk truncation below the snow line prevents the formation of giant planets (Pichardo et al. 2005). In GJ65, the snowline is located well within the stability region of the binary ~0.04 au (Ida & Lin 2005). This is a priori compatible with the formation of a Neptune type planet in this region (Table 2). The current detection thus provides an important example to study the formation mechanism of these systems, although a complete dynamical and evolutionary study of GJ65 is out of the scope of this Letter. Interestingly, GJ65AB shares strong similarities with Alpha Centauri AB from the standpoint of stability, since these two systems share almost identical mass ratio and eccentricity, which are, at first order, the sole parameters with the semi-major axis entering in $a_c$. GJ65 ($a = 5.45$ au) is scaled down by a factor 4.2 with respect to Alpha Cen AB ($a = 23.52$ au), which could form a basis for theoretical long-term stability studies of these two objects (Quarles & Lissauer 2016; Cuello & Sucerquia 2024).

The inclination of the planet's orbit can be compared to the inclination of the stars' rotation axes and that of the binary star orbit. The inclinations of the rotation axes of GJ65A and GJ65B were measured through spectro-polarimetric observations by Barnes et al. (2017), who derived inclinations of $i_A = 60 \pm 6^\circ$ and $i_B = 64 \pm 7^\circ$ ($0 < i < 90^\circ$), respectively, modulo 180°. They can therefore be either $i = 60^\circ$ or $i = 180-60 = 120^\circ$. The first scenario ($i = 60^\circ$) would indicate a ≈180° misalignment between the spin and the orbital plane of the binary, that is, a retrograde orbital configuration. This appears particularly unlikely in GJ65AB, considering that both stars are fast rotators with consistent rotation axis inclinations. We therefore favor an inclination of 120° for the stellar rotation axes, as it corresponds to a spin-orbit alignment of star and the binary. In this configuration, the inclination of the planet's orbital plane ($i_p \approx 90^\circ$; Table 2) with respect to the binary star's and the two stars' equatorial plane is $\Delta i \approx 30 \pm 10^\circ$. A similarly high relative inclination is known in several multi-planet systems between planet orbits as, for instance, $\pi$ Men (De Rosa et al. 2020; Hatzes et al. 2022), $\nu$ And (McArthur et al. 2010) or Kepler-108 (Mills & Fabrycky 2017). While the configuration is different in GJ65 with two stars A and B instead of multiple planets, the high eccentricity of their orbit ($e \approx 0.62$; Table 1) could induce the Kozai-Lidov effect, thereby raising the planet's inclination (Naoz 2016; Kozai 1962; Lidov 1962). Testing this hypothesis would benefit from follow-up observations to improve the precision of the inclination of the planetary orbit, as well as complementary measurements of the position angle and inclination of the stars' rotation axes, for example, using spectro-interferometry (Le Bouquin et al. 2009).

The immediate vicinity of GJ65 raises the question of characterizing this system with direct imaging. The age of the system was estimated to approximately ~1 Gyr from mass-luminosity isochrones by Kervella et al. (2016), based on the evolutionary models for low-mass stars in Baraffe et al. (2015). Given this system is relatively old, the intrinsic emission of the planet will be extremely low. The luminosity of the planet will be given by the reflected light contrast $c = g(\phi) \cdot A \cdot (\frac{R_p}{a})^2$ with A as the
albedo, $g(\phi)$ the phase function on the orbit, and $R_p$ the planet radius. The estimated contrast is $c = 10^{-6.3}$ at a projected separation of 90mas, assuming $g(\phi) \cdot A = 0.4$, and $R_p = 7 R_⊕$ based on the mass-radius relation of Chen & Kipping (2017). This value of contrast is at the limit of current observation capabilities, but could be a prime target for direct high-contrast observations with VLTI/GRAVITY+ and its planned extreme adaptive optics (GRAVITY+ Collaboration 2022b), as well as the upcoming ELT instruments (Davies et al. 2018; Carlotti et al. 2022; Palle et al. 2023). In addition, the knowledge of the orbital parameters of the planet would allow for an astrometry-informed pointing for GRAVITY+ and direct imaging observations from ground and space. Gaia DR4 will be useful to further constrain this system and potentially provide absolute astrometry to lift the degeneracy which of the two stars host the candidate planet. However, the astrometric performance may be more affected by flares in the visible compared to infrared observations. The monitoring of GJ65 with follow-up GRAVITY observations will be important in further constraining the orbital parameters of the candidate planet and the possible presence of additional planets. These observations demonstrate the ability to use astrometry to detect planets around low-mass stars, and could be applied to the monitoring of nearby <25 pc system to constrain the presence of giant planets down to Neptune mass. In the longer term, the extension of ground-based interferometers to micro-arcsecond precision could enable to probe Earth-mass planets in our close solar neighborhood.

**Acknowledgements.** We thank the anonymous referee for his comments which improved the content of this Letter. We are very grateful to our funding agencies (MPG, ERC, CNRS [PNCG, PNGRAM], DFG, BMBF, Paris Observatory [CS, PhyFOG], Observatoire des Sciences de l'Univers de Grenoble, and the Fundação para Ciência e Tecnologia), to ESO and the Paranal staff, and to the many scientific and technical staff members in our institutions, who helped to make GRAVITY a reality. P.K. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (projects CepBin, grant agreement 695099, and UniverseScale, grant agreement 951549). A.A. and P.G. acknowledge the financial support provided by FCT/Portugal through grants PTDC/FIS-AST/7002/2020 and UIDB/00099/2020. F.W. has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 10100471. This research has made use of the Jean-Marie Mariotti Center Aspro service. Based on observations collected at the European Southern Observatory under ESO programme(s): 099.C-2017(A), 099.C-2017(B), 099.C-2017(C), 0100.C-0597(A), 0100.C-0597(B), 0100.C-0597(C), 0100.C-0608(A), 0100.C-0608(B), 0101.C-0121(A), 0101.C-0121(B), 0101.C-0121(C), 0102.C-0205(A), 0102.C-0205(B), 0102.C-0205(C), 0102.C-0211(A), 0102.C-0211(B), 0102.C-0211(C), 0103.C-0183(A), 0103.C-0183(B), 0103.C-0183(C), 0103.C-0192(A), 0103.C-0183(A), 0103.C-0183(B), 0103.C-0183(C), 0104.C-0161(A), 0104.C-0161(B), 0104.C-0161(C), 0104.C-0161(D), 105.20KT.001, 105.20KT.002, 105.20KT.003, 106.215T.001, 106.215T.002, 106.215T.003, 109.2310.001, 109.2310.002, 109.2310.003, 109.2310.004, 110.23TC.001, 110.23TC.002, 110.23TC.003, 111.24KS.001, 111.24KS.002, 111.24KS.003, 111.24KS.004, 111.24KS.005, 111.24KS.006, 112.25GN.001, 112.25GN.002, 112.25GN.003. This research has made use of Astropy (http://www.astropy.org/), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), the Numpy library (Harris et al. 2020), the Scipy library (Virtanen et al. 2020) and the Matplotlib graphics environment (Hunter 2007).

## References

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

* Barnes, J. R., Jeffers, S. V., Haswell, C. A., et al. 2017, MNRAS, 471, 811

* Bouchy, F., Doyon, R., Artigau, É., et al. 2017, The Messenger, 169, 21

* Carlotti, A., Jocou, L., Moulin, T., et al. 2022, in Adaptive Optics Systems VIII, eds. L. Schreiber, D. Schmidt, & E. Vernet, SPIE Conf. Ser., 12185, 121855L

* Catanzarite, J., Shao, M., Tanner, A., Unwin, S., & Yu, J. 2006, PASP, 118, 1319

* Chabrier, G., & Baraffe, I. 2000, ARA&A, 38, 337

* Chen, J., & Kipping, D. 2017, ApJ, 834, 17

* Colavita, M. M., Wallace, J. K., Hines, B. E., et al. 1999, ApJ, 510, 505

* Cuello, N., & Sucerquia, M. 2024, Universe, 10, 64

* Davenport, J. R. A., Becker, A. C., Kowalski, A. F., et al. 2012, ApJ, 748, 58

* Davies, R., Alves, J., Clénet, Y., et al. 2018, in Ground-based and Airborne Instrumentation for Astronomy VII, eds. C. J. Evans, L. Simard, & H. Takami, SPIE Conf. Ser., 10702, 107021S

* De Rosa, R. J., Dawson, R., & Nielsen, E. L. 2020, A&A, 640, A73

* Delplancke, F., Derie, F., Lévêque, S., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, eds. J. D. Monnier, M. Schöller, & W. C. Danchi, 6268, 62680U

* Eisenhauer, F., Monnier, J. D., & Pfuhl, O. 2023, ARA&A, 61, 237

* Eriksson, U., & Lindegren, L. 2007, A&A, 476, 1389

* Gaia Collaboration (Arenou, F., et al.) 2023, A&A, 674, A34

* Gillessen, S., Lippa, M., Eisenhauer, F., et al. 2012, in Optical and Infrared Interferometry III, eds. F. Delplancke, J. K. Rajagopal, & F. Malbet, SPIE Conf. Ser., 8445, 844510

* GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

* GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 618, L10

* GRAVITY Collaboration (Lacour, S., et al.) 2019, A&A, 623, L11

* GRAVITY Collaboration (Abuter, R., et al.) 2022a, A&A, 657, A82

* GRAVITY+ Collaboration (Abuter, R., et al.) 2022b, The Messenger, 189, 17

* Halbwachs, J.-L., Pourbaix, D., Arenou, F., et al. 2023, A&A, 674, A9

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Hatzes, A. P., Gandolfi, D., Korth, J., et al. 2022, AJ, 163, 223

* Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621

* Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

* Ida, S., & Lin, D. N. C. 2005, ApJ, 626, 1045

* Joy, A. H., & Humason, M. L. 1949, PASP, 61, 133

* Kervella, P., Mérand, A., Ledoux, C., Demory, B. O., & Le Bouquin, J. B. 2016, A&A, 593, A127

* Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72

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

* Kozai, Y. 1962, AJ, 67, 591

* Lacour, S., Eisenhauer, F., Gillessen, S., et al. 2014, A&A, 567, A75

* Lagrange, A. M., Meunier, N., Desort, M., & Malbet, F. 2011, A&A, 528, L9

* Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, in Optical and Infrared Interferometry IV, eds. J. K. Rajagopal, M. J. Creech-Eakman, & F. Malbet, SPIE Conf. Ser., 9146, 91462D

* Le Bouquin, J. B., Absil, O., Benisty, M., et al. 2009, A&A, 498, L41

* Lidov, M. L. 1962, Planet. Space Sci., 9, 719

* Lippa, M., Gillessen, S., Blind, N., et al. 2016, in Optical and Infrared Interferometry and Imaging V, eds. F. Malbet, M. J. Creech-Eakman, & P. G. Tuthill, SPIE Conf. Ser., 9907, 990722

* Luyten, W. J. 1949, ApJ, 109, 532

* Marzari, F., & Thebault, P. 2019, Galaxies, 7, 84

* McArthur, B. E., Benedict, G. F., Barnes, R., et al. 2010, ApJ, 715, 1203

* Meunier, N. 2021, arXiv e-prints [arXiv:2104.06072]

* Mills, S. M., & Fabrycky, D. C. 2017, AJ, 153, 45

* Muterspaugh, M. W., Lane, B. F., Kulkarni, S. R., et al. 2010, AJ, 140, 1657

* Naoz, S. 2016, ARA&A, 54, 441

* Palle, E., Biazzo, K., Bolmont, E., et al. 2023, Exp. Astron., submitted [arXiv:2311.17075]

* Pfuhl, O., Haug, M., Eisenhauer, F., et al. 2014, in Optical and Infrared Interferometry IV, eds. J. K. Rajagopal, M. J. Creech-Eakman, & F. Malbet, SPIE Conf. Ser., 9146, 914623

* Pichardo, B., Sparke, L. S., & Aguilar, L. A. 2005, MNRAS, 359, 521

* Quarles, B., & Lissauer, J. J. 2016, AJ, 151, 111

* Shao, M., & Colavita, M. M. 1992, A&A, 262, 353

* Unwin, S. C., Shao, M., Tanner, A. M., et al. 2008, PASP, 120, 38

* Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261

* Winters, J. G., Henry, T. J., Jao, W.-C., et al. 2019, AJ, 157, 216

* Waillez, J., & Lacour, S. 2013, ApJ, 764, 109

* Waillez, J., Akeson, R., Colavita, M., et al. 2010, in Optical and Infrared Interferometry II, eds. W. C. Danchi, F. Delplancke, & J. K. Rajagopal, SPIE Conf. Ser., 7734, 773412

* Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

<footnote n="1">1 Max Planck Institute for Extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany</footnote>

<footnote n="2">2 LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 Place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="3">3 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>
⁴ European Southern Observatory, Karl-Schwarzschild-Straße 2,
85748 Garching, Germany

⁵ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

⁶ 1st Institute of Physics, University of Cologne, Zülpicher Straße 77,
50937 Cologne, Germany

⁷ CENTRA – Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal

⁸ Universidade de Lisboa – Faculdade de Ciências, Campo Grande,
1749-016 Lisboa, Portugal

⁹ European Southern Observatory, Casilla 19001, Santiago 19,
Chile

¹⁰ Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

¹¹ Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112,
21029 Hamburg, Germany

¹² Departments of Physics & Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA

¹³ Max Planck Institute for Astrophysics, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany

¹⁴ Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany

¹⁵ ORIGINS Excellence Cluster, Boltzmannstraße 2, 85748 Garching, Germany

¹⁶ Department of Physics, Technical University of Munich, 85748 Garching, Germany

¹⁷ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, D02 XF86 Dublin, Ireland

¹⁸ School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

¹⁹ Institute of Astronomy, University of Cambridge, Madingley Rd., Cambridge CB3 0HA, UK
# Appendix A: Swap observations

The goal of the swap sequence is to measure the phases of Eq. 1 with two opposite signs for the field component, while keeping constant the zero-point of the metrology. This is done using the internal K-mirror of GRAVITY fiber coupler (Pfuhl et al. 2014), which allows us to rotate the field before the injection of light in the fibers, given the separation of the stars and the fact that both components are equally bright and can be equally used as a fringe-tracking star.

The two phases which are measured for each step of the swap sequence are:

$$ \begin{cases} \phi_1 = \frac{2\pi}{\lambda}(\vec{s}_B - \vec{s}_A) \cdot \vec{B}_{NAB,1} + \phi_0(\lambda) \\ \phi_2 = \frac{2\pi}{\lambda}(\vec{s}_A - \vec{s}_B) \cdot \vec{B}_{NAB,2} + \phi_0(\lambda) \end{cases} \qquad (A.1) $$

For each baseline, the metrology zero point can finally be extracted from the combination of Eq. A.1 pointings:

$$ \begin{aligned} \phi_0(\lambda) &= \frac{1}{2} \arg \left[ \tilde{V}_{FT,1} e^{-\frac{2\pi}{\lambda}(\vec{s}_B - \vec{s}_A)\vec{B}_{NAB,1}} \cdot \tilde{V}_{FT,2} e^{-\frac{2\pi}{\lambda}(\vec{s}_A - \vec{s}_B)\vec{B}_{NAB,2}} \right] \\ &= \frac{1}{2} \arg \left[ \tilde{V}_{FT,1} \cdot \tilde{V}_{FT,2} \right] - \frac{\pi}{\lambda} (\vec{s}_B - \vec{s}_A) \cdot (\vec{B}_{NAB,1} - \vec{B}_{NAB,2}) \end{aligned} \qquad (A.2) $$

with $\tilde{V}_{FT,k}$ the complex visibility amplitude measured by the fringe-tracker. The metrology zero-point is then included in Eq. A.1, which provides the relative separation $\vec{s}_B - \vec{s}_A$ of components A and B.

# Appendix B: Orbit fitting

The stellar orbit and the planetary orbits are both defined using the Thiele-Innes elements, related to the four orbital elements $(a, e, \omega, \Omega)$. We followed the convention in Gaia Collaboration (2023), Halbwachs et al. (2023):

$$ \begin{cases} A = \cos \Omega \cos \omega - \sin \Omega \sin \omega \cos i \\ B = \sin \Omega \cos \omega + \cos \Omega \sin \omega \cos i \\ F = -\cos \Omega \sin \omega - \sin \Omega \cos \omega \cos i \\ G = -\sin \Omega \sin \omega + \cos \Omega \cos \omega \cos i \end{cases} \qquad (B.1) $$

The projected separation of the secondary with respect to the primary is then defined as:

$$ \begin{cases} \Delta\delta = a(AX + FY) \\ \Delta\alpha^* = a(BX + GY) \end{cases} \qquad (B.2) $$

where X and Y are the elliptical rectangular coordinates:

$$ \begin{cases} X = \cos E - e \\ Y = (\sin E) \sqrt{1 - e^2} \end{cases} \qquad (B.3) $$

with E as the eccentric anomaly, obtained by solving the Kepler equation.

<figure><img src="image_5.png" /><figcaption>Fig. B.1. Period search of GJ65 residuals, with all orbital parameters free for each (fixed) period and including covariances, shown with inverse reduced chi-square in y-axis. The best period, $p = 156$ d, coincides with the peak of the Lomb-Scargle periodogram.</figcaption></figure>

# Appendix C: Stability domain

We use the stability criterion of test particles derived in Holman & Wiegert (1999) (Eq. 1), which we reproduce here:

$$ a_c = [0.464 - 0.380\mu - 0.631e + 0.586\mu e + 0.150e^2 - 0.198\mu e^2]a_b \qquad (C.1) $$

with $a_c$, as the critical semi-major axis defining the outer edge of the stability region, $\mu = m_B/m_A$ the mass ratio of the binary, $e$ the eccentricity of the binary, and $a_b$ its semi-major axis. The Neptune-mass candidate as well as the stability region of the system are shown in Figure C.1.

<figure><img src="image_6.png" /><figcaption>Fig. C.1. Parameter space probed by GRAVITY observations. The stability region is computed from Holman & Wiegert (1999).</figcaption></figure>
Appendix D: System architecture

<figure><img src="image_7.png" /><figcaption>**Fig. D.1.** Planet candidate and barycentric ellipsis of the A and B stars in GJ65. The gray-shaded disk corresponds to the region of stability, shown in the orbital plane of the binary. The blue orbit corresponds to the planet in either star A and star B hypothesis. The inset shows the on-sky orbit around GJ65 B and the location of the ice-line (green) in the orbital plane of the planetary system.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. D.2. Relative angular separation of the stars A and B in GJ65.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. D.3. Relative radial velocity of stars A and B in GJ65.</figcaption></figure>