Chasing extreme planetary architectures

G. Chauvin, M. Videla, H. Beust, R. Mendez, A. Correia, S. Lacour, A.

Tokovinin, J. Hagelberg, F. Bouchy, I. Boisse, et al.

► To cite this version:

G. Chauvin, M. Videla, H. Beust, R. Mendez, A. Correia, et al. Chasing extreme planetary architectures. Astronomy & Astrophysics - A&A, 2023, 675, pp.A114. 10.1051/0004-6361/202244502. hal-04309335

HAL Id: hal-04309335

https://cnrs.hal.science/hal-04309335v1

Submitted on 1 Dec 2023

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>
# Chasing extreme planetary architectures

## I. HD 196885 Ab, a super-Jupiter dancing with two stars?

G. Chauvin<sup>1,2</sup>®, M. Videla<sup>3,1</sup>®, H. Beust<sup>4</sup>, R. Mendez<sup>5</sup>®, A. C. M. Correia<sup>6,7</sup>®, S. Lacour<sup>8</sup>®, A. Tokovinin<sup>9</sup>®, J. Hagelberg<sup>10</sup>®, F. Bouchy<sup>10</sup>, I. Boisse<sup>11</sup>, C. Villegas<sup>3,1</sup>, M. Bonavita<sup>12</sup>®, S. Desidera<sup>13</sup>®, V. Faramaz<sup>14</sup>, T. Forveille<sup>4</sup>®, A. Gallenne<sup>15,1</sup>®, X. Haubois<sup>16</sup>®, J. S. Jenkins<sup>17,18</sup>, P. Kervella<sup>8</sup>®, A.-M. Lagrange<sup>8</sup>®, C. Melo<sup>16</sup>, P. Thebault<sup>8</sup>, S. Udry<sup>10</sup>, and D. Segransan<sup>10</sup>

(Affiliations can be found after the references)

Received 15 July 2022 / Accepted 22 September 2022

## ABSTRACT

**Context.** Planet(s) in binaries are unique architectures for testing predictions of planetary formation and evolution theories in very hostile environments. Their presence in such a highly perturbed region poses a clear challenge from the early phase of planetesimals accretion to the dynamical evolution and stability through a very long period of time (several Gyrs in some case).

**Aims.** The combination of radial velocity, speckle interferometry, high-contrast imaging and high-precision astrometry with interferometry, offers a unique and unprecedented set of observable to push the exploration of the extreme planetary system HD 196885. By dissecting the physical and orbital properties of each component, we aim at shedding light on its global architecture and stability.

**Methods.** We used the IRDIS dual-band imager of SPHERE at VLT, and the speckle interferometric camera HRCAM of SOAR, to acquire high-angular resolution images of HD 196885 AB between 2015 and 2020. Radial velocity observations started in 1983 with CORAVEL at OHP have been extended with a continuous monitoring with CORALIE at La Silla, and ELODIE and SOPHIE at OHP over almost 40 yr extending the radial velocity measurements HD 196885 A and resolving both the binary companion and the inner giant planet HD 196885 Ab. Finally, we took advantage of the exquisite astrometric precision of the dual-field mode of VLTI/GRAVITY (down to 30 µas) to monitor the relative position of HD 196885 A and B to search for the 3.6 yr astrometric wobble of the circumpri-mary planet Ab imprinted on the binary separation.

**Results.** Our observations enable to accurately constrain the orbital properties of the binary HD 196885 AB, seen on an inclined and retrograde orbit ($i_{AB} = 120.43$ deg) with a semi-major axis of 19.78 au, and an eccentricity of 0.417. The GRAVITY measurements con-firm for the first time the nature of the inner planet HD 196885 Ab by rejecting all families of pole-on solutions in the stellar or brown dwarf masses. The most favored island of solutions is associated with a Jupiter-like planet ($M_{Ab} = 3.39 M_{Jup}$), with moderate eccentricity ($e_{AaAb} = 0.44$), and inclination close to 143.04 deg. This results points toward a significant mutual inclination ($\Phi = 24.36$ deg) between the orbital planes (relative to the star) of the binary companion B and the planet Ab. Our dynamical simulations indicate that the system is dynamically stable over time. Eccentricity and mutual inclination variations could be expected for moderate von Zipele Kozai Lidov cycles that may affect the inner planet.

**Key words.** planet-star interactions – binaries close – planets and satellites: dynamical evolution and stability – instrumentation: high angular resolution – instrumentation: adaptive optics – instrumentation: spectrographs

## 1. Introduction

Even though binaries were initially banished from exoplanet hunting campaigns, more than 200 planets<fnref n="1" /> have been so far discovered in multiple systems orbiting one component or both components, referring to the S-type and P-type configuration respectively (see classification by Dvorak 1982). These systems have been initially discovered by accident like Gliese 86 Ab (Queloz et al. 2000) and γ Cep Ab (Hatzes et al. 2003), but more recently, dedicated surveys in radial velocity (RV), transit and imaging have been targeting them (e.g. Fontanive & Bardalez Gagliuffi 2021). Given that the formation of stars in multiple systems is a frequent byproduct of stellar formation, a natural question is to understand how the presence of a stellar companion can affect the planetary formation process. Planetary formation theories have been restricted for years to the case

of a single star environment to understand the formation of our own Solar System (Mordasini et al. 2008), and more recently investigated for binary stars (e.g. Silsbee & Rafikov 2021).

For the most frequent dynamical configuration observed for planet(s) in binaries, the S-type circumprimary one with a planet orbiting one component of the binary (generally the most massive one), models have predicted that the presence of a very close binary companion can hinder planetary formation around a primary star in several ways. The most problematic issue is probably not the disk truncation (see Jang-Condell 2015), but rather the dynamical excitation that could prevent mutual planetesimal accretion for objects in the 100 m-to-100 km size range (Thebault et al. 2006; Thebault 2011; Paardekooper et al. 2008; Fragner et al. 2011), hence obstructing the formation of a planet by core accretion, or ejecting the planet in unstable systems (e.g. Thebault & Haghighipour 2015). Dedicated observing campaigns suggest a transition at typical separation within 50–300 au indicating that short or intermediate-period binaries have statistically less chance to host planets or brown dwarf companions,

<footnote n="1">¹ A total of 217 planets in binaries distributed over 154 systems are known on July, 2021 (https://www.univie.ac.at/adg/schwarz/multiple.html).</footnote>
compared to wide binaries that on the contrary have no influ-
ence on the architectures of planetary systems (Ziegler et al.
2021). This effect seems to be corroborated by the study of young
stars for which short-separation (≤100 au) binaries have a lower
probability of hosting circumstellar dust in the innermost few au
around each star, therefore with a depleted reservoir of solids for
the formation of planets by core accretion (Duchêne 2010).

In the last decade, more intensive and dedicated surveys have
been dedicated to search for S-type and P-type<fnref n="2" /> exoplanetary sys-
tems using (and often combining) RV, transit, high-contrast and
direct imaging observations. RV surveys, able to tackle multiple
stellar signals at the precision required for detecting planets, were
specifically carried out (Konacki 2005; Eggenberger & Udry
2007; Toyota et al. 2009; Desidera et al. 2010). Zucker & Mazeh
(2002) initially showed that the planet properties in binaries were
different to planets around single stars. Eggenberger et al. (2004)
confirmed that the most massive short period planets were found
in binary systems and that planets with short (P ≤ 40 days) peri-
ods in binaries were likely to have low eccentricities. Tamuz et al.
(2008) finally found that extremely high eccentric planets were
all in binaries.

In parallel to Doppler surveys, high-contrast imaging
campaigns have been able to access intermediate-separation
(20–100 au) binaries where the companion influence is the
most expected. Several deep imaging studies have revealed a
number of additional close companions to known exoplanetary
hosts (Patience et al. 2002; Luhman & Jayawardhana 2002;
Chauvin et al. 2006; Eggenberger & Udry 2007; Law et al. 2014;
Mugrauer & Ginski 2015; Baranec et al. 2016; Ginski et al.
2021). Eggenberger et al. (2007) carried out a dedicated survey to
explore the multiplicity in two samples of stars with and without
planets detected in RV in order to test the impact of duplicity on
the giant planet occurrence. They found a lower corrected binary
fraction for the planet-host sub-sample, particularly for physical
separations shorter than 100 au (Eggenberger et al. 2008).

For longer period (≥100 au) binaries, catalogs compiled
by Raghavan et al. (2006), Desidera & Barbieri (2007), and
Bonavita & Desidera (2007) further updated by Bonavita &
Desidera (2020) studied the effect of duplicity, and did not find
any significant discrepancies in the planet frequency between
binary and single stellar systems. More recently, the use of
the Gaia offered the opportunity of a systematic search for
(sub)-stellar companions to exoplanet host (Kervella et al.
2019; Mugrauer 2019; Michel & Mugrauer 2021). Fontanive &
Bardalez Gagliuffi (2021) recently derived an overall raw mul-
tiplicity rate of 23.2±1.6% for hosts to exoplanets and brown
dwarf companions, considering a volume-limited sample up to
200 pc cross-matched with the Gaia Data Release 2 catalog
(Gaia Collaboration 2016a, 2018).

A complementary approach has been to target individual
systems, and particularly the extreme ones that provide theo-
reticians with an ideal test bed for the theories of planetary
formation, evolution and dynamical stability in a very hostile
environment. Today, there are about 16 extreme exoplanetary
systems (binaries with small semi-major axes of less than about
20 au hosting a planet; see Fig. 1). One emblematic case is γ Cep
for which Benedict et al. (2018) managed to derive the mass and
the orbital properties of the circumprimary planet using accu-
rate astrometric measurements with HST. Orbiting at 1.94 au in
a highly mutually inclined orbit (Φ = 70 ± 13 deg), the formation
and survival of this 9.4<sup>+0.7</sup><sub>-1.1</sub> M<sub>Jup</sub> exoplanet remains unresolved

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Overview of the binary star systems. The sizes of the cir-
cles representing the primary and secondary star indicate their masses
according to the scale reported on the *bottom-right*. The planet masses
are not to scale. The green zone marks the region of orbital stability for
a coplanar configuration in the context of an S-type configuration.</figcaption></figure>

Bazsó et al. (2017). Another interesting system is HD 196885 that Chauvin et al. (2006) discovered in the course of a VLT/CFHT deep imaging survey of exoplanet hosts. In this paper, we present the results of the long-term monitoring of this system combining RV, direct imaging, speckle and dual-field interferometry. The physical properties of the HD 196885 exoplanetary system are reported in Sect. 2. The observations and data reduction are presented in Sect. 3, together with the results in Sect. 4. In Sect. 5, we finally discuss their implication in terms of system architecture, formation and stability.

## 2. The HD 196885 exoplanetary system

HD 196885 is an F8V (V = 6.385 mag, H = 5.19 mag) star
located at 34.00 ± 0.04 pc (Gaia EDR3, Gaia Collaboration
2021). Based on CORALIE spectra, Sousa et al. (2006) derived
a spectroscopic temperature, surface gravity, and metallicity
of T<sub>eff</sub> = 6340 ± 39 K, log(g) = 4.46 ± 0.02 and [Fe/H] =
0.29 ± 0.05, respectively. The projected stellar rotational velocity
v sin(i) was estimated to 7.3 ± 1.5 km s<sup>-1</sup> from ELODIE spectra
(Correia et al. 2008). These results were supported by indepen-
dent Lick observations reported by Fischer et al. (2009). The
chromospheric activity level of log(R'<sub>HK</sub>) = −5.01 (Wright et al.
2004) is relatively low when placed into context with the activ-
ity distribution of the solar neighborhood (Jenkins et al. 2011).
Bolometric luminosity correction and evolutionary model pre-
dictions lead to an estimate of the luminosity and the mass of
2.4 L<sub>☉</sub> and 1.3 M<sub>☉</sub>, respectively. The corresponding stellar age
derived from evolutionary tracks and from the activity level
varies between 1.5 and 3.5 Gyr (see Correia et al. 2008; Fischer
et al. 2009). Finally, we looked at wider physical companions
within 3 arcmin from HD 196895 in Gaia EDR3 and found none.

<footnote n="2">² P-type orbits occur when the planet orbits both binary components, and generally referred as circumbinary planet.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** SPHERE/IRDIS observations of HD 196885 AB obtained in September, 2017. The early-M1V dwarf companion HD 196885 B is resolved at *K1*-band at a separation of 604 mas and a position angle of 45.7 deg.</figcaption></figure>

The companion BUP 1908 (192 arcseconds separation, position angle 8°, $V = 10.4$ mag, TYC 1096-637-1) listed in the double star catalog, CCDM, is optical according to Gaia. As expected, HD 196885 shows a huge signature of binarity (with an S/N of 221) in the Hipparcos-Gaia proper motion anomaly catalog (Kervella et al. 2022), but less constraining that the set of observations presented below.

In 2004, an RV variation indicating a planetary candidate was measured at Lick observatory and adjusted with a preliminary orbital period of $P = 0.95$ yr. The result was temporarily reported on the California Planet Search Exoplanet Web site, but finally withdrawn owing a significant residual drift in the orbital solution. Nevertheless, this star was included in a VLT/CFHT deep-imaging survey of stars hosting planets detected by RV observations (Chauvin et al. 2006). This led them to the discovery of a close (~20 au) stellar binary companion to HD 196885, motivating dedicated follow-up observations with Adaptive-Optics (AO) imaging. The recent image of the system obtained by SPHERE at VLT is shown in Fig. 2. Following the CFHT discovery, astrometric monitoring and long-slit spectroscopic observations at VLT with NaCo enabled to resolve the orbital motion of the stellar companion and to confirm its low stellar-mass nature as an $M_{1\pm1V}$ dwarf companion located at only 0.7 as (23 au in projected physical separation) and likely to be responsible for the trend seen in the Lick RV residuals (Chauvin et al. 2007). Using a double-Keplerian model for the binary star and the planet to adjust their ELODIE, CORALIE, and CORAVEL observations spread over 14 yr, Correia et al. (2008) confirmed the existence of the planet HD 196885 Ab orbiting the primary. They derived a first range of orbital solutions with a minimum mass of $M_{Ab} \sin(i_{Ab}) = 2.96 M_{Jup}$, a period of $P_{Ab} = 3.69 \pm 0.03$ yr, and an eccentricity of $e_{AbAb} = 0.462 \pm 0.026$. Moreover, they found additional constraints for the binary companion HD 196885 B with a period of $P_{AB} \ge 40$ yr, a semi-major axis $a_{AB} \ge 14$ au, and a minimum mass of $M_B \sin(i_{AB}) \ge 0.28 M_{\oplus}$. Based on Lick observations,

Fischer et al. (2009) derived consistent results for both the inner planet and the binary companion.

In 2011, Chauvin et al. (2011) reported a 4 yr astrometric monitoring of the separation between HD 196885 A and B using NaCo at VLT. Combined with RV observations, these new measurements enabled them to refine the constraints on the orbital properties of the stellar companion, and to derive an inclination of $i_{AB} = 116.8 \pm 0.7$ deg, an eccentricity of $e_{AB} = 0.42 \pm 0.03$, a semi-major axis of $a_{AB} = 21.00 \pm 0.86$ au, and a true mass of of $M_B = 0.45 \pm 0.01 M_{\oplus}$ consistent with the spectral type M1±1V derived by Chauvin et al. (2007). The key physical properties of the extreme planetary system HD 196885, the stellar host (Aa), the binary companion (B) and the circumbinary planet (Ab) based on the literature can be found in Table 1.

# 3. Observations and data reduction

## 3.1. Radial velocities

In an attempt to revisit the physical and orbital properties of the HD 196885 system, we considered all data from past radial velocity surveys already published in Chauvin et al. (2011) for the primary host HD 196885 A: CORAVEL (9 RV measurements from June 1982 to August 1997), ELODIE (69 measurements from June 1997 to August 2006), CORALIE (33 measurements from April 1999 to November 2002), and Lick (75 measurements from 1998 to 2008). To this, we have added new SOPHIE observations (52 measurements from 2007 to July 2017) to gather a total of 238 radial velocity measurements spanning almost 40 yr of observation (see Correia et al. 2008; Fischer et al. 2009; Chauvin et al. 2011, and this work). SOPHIE data were acquired in High-resolution mode as part of the program dedicated to the follow-up of ELODIE candidates. Data from 2007 to June 2011 were reduced thanks to the SOPHIE pipeline (Bouchy et al. 2009). In June 2011, an update of the instrument was done. It consists mainly in the addition of octagonal fibers to light guiding from the telescope to the spectrograph (Bouchy et al. 2013). Data after this date were then reduced as in Heidari et al. (2022).

## 3.2. SPHERE adaptive-optics imaging

Following previous astrometric observations at VLT obtained with the NaCo instrument, the HD 196885 system has been monitored during the SpHere Infrared survey for Exoplanets (SHINE, Chauvin et al. 2017; Desidera et al. 2021; Langlois et al. 2021) at 2 different epochs on 5 June 2015 and 30 September 2017 (see Table 2) using the VLT/SPHERE high-contrast instrument (Beuzit et al. 2019). The observations were obtained with the modes IRDIFS-EXT that combine simultaneously the IRDIS (Dohlen et al. 2008) and IFS instruments (Claudi et al. 2008), although IRDIS observations were mainly used in the current study. The IRDIFS-EXT mode combines IRDIS in dual-band imaging (DBI, Vigan et al. 2010) mode with the K1K2 filter doublet $\lambda_{K1} = 2.103 \pm 0.102$ µm, $\lambda_{K2} = 2.255 \pm 0.109$ µm, and IFS in the YJH (0.97–1.66 µm) setting. Each observing sequence was performed with the field-tracking mode. The detail about the observations are reported in Table 2. Only the IRDIS observations at *K1*-band (with a lower background noise) were used to extract the relative astrometry of HD 196885 B relative to A.

## 3.3. SOAR speckle interferometry

Complementary astrometric measurements of HD 196885 have been obtained with the High-Resolution Camera (HRCAM; Tokovinin et al. 2021) at the SOAR 4.1-m telescope
<table><caption>**Table 1.** Physical properties of the extreme planetary system HD 196885, the stellar host (Aa), the binary companion (B) and the circumprimary planet (Ab) based on the latest studies from Correia et al. (2008); Fischer et al. (2009); Chauvin et al. (2011); Gaia Collaboration (2021).</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td colspan="2">Primary star: HD 196885 Aa</td></tr><tr><td>Spectral type</td><td>F8V</td></tr><tr><td>G (mag)</td><td>6.276973 ± 0.002767</td></tr><tr><td>V (mag)</td><td>6.385 ± 0.010</td></tr><tr><td>H (mag)</td><td>5.190 ± 0.023</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>6340 ± 39</td></tr><tr><td>log(g) (dex)</td><td>4.46 ± 0.02</td></tr><tr><td>Fe/H (dex)</td><td>0.29 ± 0.05</td></tr><tr><td>Distance (pc)</td><td>34.00 ± 0.04</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>71.915 ± 0.032</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>89.318 ± 0.026</td></tr><tr><td>Age (Gyr)</td><td>1.5–3.5</td></tr><tr><td>M<sub>Aa</sub> (M<sub>☉</sub>)</td><td>1.3 ± 0.1</td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td>2.4 ± 0.1</td></tr><tr><td colspan="2">Exoplanet: HD 196885 Ab</td></tr><tr><td>M<sub>Ab</sub>sin(i<sub>AaAb</sub>) (M<sub>Jup</sub>)</td><td>2.98 ± 0.05</td></tr><tr><td>P<sub>AaAb</sub> (yr)</td><td>3.63 ± 0.01</td></tr><tr><td>e<sub>AaAb</sub></td><td>0.48 ± 0.02</td></tr><tr><td>ω<sub>AaAb</sub> (deg)</td><td>93.2 ± 3.0</td></tr><tr><td>i<sub>P_AaAb</sub></td><td>2002.85 ± 0.02</td></tr><tr><td>a<sub>AaAb</sub> (au)</td><td>2.6 ± 0.1</td></tr><tr><td colspan="2">Binary companion: HD 196885 B</td></tr><tr><td>Spectral type</td><td>M1 ± 1</td></tr><tr><td>K<sub>B</sub> (mag)</td><td>6.398 ± 250</td></tr><tr><td>M<sub>B</sub> (M<sub>☉</sub>)</td><td>0.45 ± 0.01</td></tr><tr><td>P<sub>AB</sub> (yr)</td><td>72.06 ± 4.59</td></tr><tr><td>e<sub>AB</sub></td><td>0.42 ± 0.03</td></tr><tr><td>ω<sub>AB</sub> (deg)</td><td>-118.1 ± 3.1</td></tr><tr><td>Ω<sub>AB</sub> (deg)</td><td>79.8 ± 0.1</td></tr><tr><td>i<sub>AB</sub> (deg)</td><td>116.8 ± 0.7</td></tr><tr><td>i<sub>P,AB</sub></td><td>1985.59 ± 0.39</td></tr><tr><td>a<sub>AB</sub> (au)</td><td>21.00 ± 0.86</td></tr></tbody></table>

**Notes.** We consider below that the system is composed of the inner binary subsystem AaAb (composed by the primary star Aa and the exoplanet Ab), and the outer binary subsystem AB (composed by the center of mass of the inner orbit A and the companion star B).

Tokovinin & Cantarutti (2008)<fnref n="3" />. While regular observations are made in the Stroemgren *y*, Cousins *I*, or Hα filters, all observations of HD 196885 were done with the *I* filter with a resolution of 36 mas. The data have been processed by a custom IDL pipeline developed by A. Tokovinin to extract the binary parameters (separation, position angle, and magnitude difference) by fitting the secondary peaks of the auto-correlation functions or the fringes from the power spectrum of the observations obtained by Fourier transform as explained in Tokovinin et al. (2010). HRCAM routinely delivers a precision of 1–3 mas in angular separation for objects brighter than V ~ 12. On every observing night “calibration binaries” are included, these are binaries with very well known orbits (grade 5 on the USNO orbit catalogue), from which the measurements are calibrated, leading to

systematic errors less than 0.1 deg in position angle, and better than 0.2% in scale, i.e., smaller than the internal precision. The final precision of these measurements depends on a number of factors, but in this paper we adopt a uniform uncertainty of 3 mas, which is representative of all the measurements (Tokovinin 2018). The Observing log and the astrometric results are reported in Tables 2 and 3.

### 3.4. VLTI/GRAVITY dual-field interferometry

HD 196885 was also monitored in the context of monitoring program with GRAVITY (Gravity Collaboration 2017) and the UTs in MED resolution at K-band. Over 2 yr, we have obtained a total of four observations on 15 August 2019, 29 May 2021, 25 August 2021, and 25 September 2021. The weather conditions ranged from good to very good, and the observing log of the observations is also reported in Table 2.

These observations were obtained with the astrometric dual-field mode in off-axis: in the first exposure, the fringe tracker single-mode fiber is centered on the primary star HD 196885 A while the scientific single-mode fiber is centered on the stellar companion HD 196885 B. In a second exposure, the two fibers are swapped, to calibrate the metrology offset. This observation mode is different from the dual-field mode in on-axis, where the planet is bright enough to be directly observed (Gravity Collaboration 2019). Here, the GRAVITY interferometer only observed the A and B components. The observations were pre-processed with the Public Release 1.5.0 (1 July 2021) of the ESO Gravity pipeline (Lapeyrere et al. 2014). The reduced data are then processed to extract astrometry using the GRAVITY python tools offered to the community<fnref n="4" />. The astrometric values are presented in Table 3.

## 4. Resolving the astrometric wobble due to the giant planet HD 196885 Ab

The main objective of our GRAVITY program was to measure the astrometric wobble in the separation of a binary system related to presence of the giant planet, and in such a circumstance be able to derive the orbital properties of the binary companion and the giant planets, but also to confirm the planetary nature of HD 196885 Ab. To date, such observation have only been performed in the context of the γ Cep system by Benedict et al. (2018). Combining all the available measurements obtained from radial velocity, adaptive-optics imaging, speckle interferometry, and dual-field interferometry, we have revisited the system properties using a Bayes-based orbital elements estimation inspired in a previous study by our team to triple hierarchical stellar systems (Villegas et al. 2021).

### 4.1. Code description and set up

The orbital parameters estimation code has been developed in a Bayesian Markov-chain Monte Carlo-based framework, in the context of the so-called hierarchical approximation to combine astrometric with spectroscopic measurements for multiple systems. The current implementation of this code<fnref n="5" /> is an adaptation

<footnote n="3">³ For up-to-date details of the instrument see http://www.ctio.noirlab.edu/~atokovin/speckle/</footnote>

<footnote n="4">⁴ https://version-lesia.obspm.fr:/repos/DRS_gravity/gravi_tools3</footnote>

<footnote n="5">⁵ The original software is available on GitHub: BinaryStars codebase: https://github.com/mvidela31/BinaryStars under a 3-Clause BSD License. Its evolution used in this study will soon be released.</footnote>
<table><caption>**Table 2.** Observation log of the new SPHERE, SOAR, and GRAVITY observations.</caption><thead><tr><th>UT-Date</th><th>MJD</th><th>Telescope</th><th>Instr.</th><th>Filter</th><th>Mode</th><th>Nb Exp.</th><th>Nb Integration</th><th>Exp. Time<br>(s)</th><th>ω<br>(")</th><th>τ<sub>0</sub><br>(ms)</th></tr></thead><tbody><tr><td>27-09-2000</td><td>51826.186</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>3</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>27-09-2002</td><td>52459.690</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>3</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>05-06-2015</td><td>57178.378</td><td>VLT/UT3</td><td>IRDIS</td><td>K1</td><td>DBI</td><td>21</td><td>35</td><td>2</td><td>1.75</td><td>1.2</td></tr><tr><td>27-09-2015</td><td>57291.493</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>3</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>30-09-2017</td><td>58026.077</td><td>VLT/UT3</td><td>IRDIS</td><td>K1</td><td>DBI</td><td>2</td><td>30</td><td>2</td><td>0.72</td><td>3.4</td></tr><tr><td>27-05-2018</td><td>58265.219</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>2</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>16-07-2019</td><td>58680.881</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>2</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>15-08-2019</td><td>58710.188</td><td>VLTI/UTs</td><td>GRAVITY</td><td>K</td><td>MED</td><td>3</td><td>20</td><td>10</td><td>0.95</td><td>2.9</td></tr><tr><td>29-05-2021</td><td>59363.331</td><td>VLTI/UTs</td><td>GRAVITY</td><td>K</td><td>MED</td><td>3</td><td>12</td><td>30</td><td>0.64</td><td>5.0</td></tr><tr><td>24-08-2021</td><td>59450.195</td><td>VLTI/UTs</td><td>GRAVITY</td><td>K</td><td>MED</td><td>3</td><td>12</td><td>30</td><td>0.90</td><td>2.6</td></tr><tr><td>02-10-2021</td><td>59489.151</td><td>SOAR</td><td>HRCAM</td><td>I</td><td>SI</td><td>3</td><td>400</td><td>0.024</td><td></td><td></td></tr><tr><td>20-10-2021</td><td>59507.082</td><td>VLTI/UTs</td><td>GRAVITY</td><td>K</td><td>MED</td><td>3</td><td>12</td><td>30</td><td>1.00</td><td>2.8</td></tr></tbody></table>

**Notes.** The exposure time is given by Nexp, the number of exposures, multiplied by NDIT and DIT, which are the number and duration of integrations per exposure (in the case of HRCAM in ms). ω and τ₀ are the seeing and the coherence time respectively. For the modes, DBI refers to Dual-Band Imaging, SI to Speckle Interferometry and MED for medium spectral resolution interferometry.

of the Bayesian inference algorithm for binary stars developed by Videla et al. (2022), which is based on the No-U-Turn Sampler (NUTS) algorithm (Hoffman et al. 2014). The NUTS algorithm is an MCMC method that avoids the random-walk behavior and the sensitivity to correlated parameters of other MCMC methods commonly used in astronomy (such as the Metropolis-Hastings within Gibbs sampler Ford 2005, the Parallel Tempering sampler Gregory 2005, 2011, the Affine Invariant MCMC Ensemble sampler Hou et al. 2012, and the Differential Evolution Markov Chain sampler Nelson et al. 2013) by incorporating first-order gradient information of the parameters space to guide the sampling steps with an adaptive criterion for determining their lengths. The NUTS algorithm qualities allow for performing a quick and efficient estimation of the target distribution, specially in the case of the high-dimensionality problem in hand. We estimate the posterior distribution of the joint orbital parameter space of a triple hierarchical system that we adapted to explore the planetary architecture in S-type binaries. The parameter space consists of the orbital parameters of the inner binary subsystem AaAb (composed by the primary star Aa and the exoplanet Ab), and the orbital parameters of the outer binary subsystem AB (composed by the center of mass of the inner orbit A and the companion star B), assuming of course the same value of the parallax π for both subsystems. Therefore, the parameter space is defined by the vector $\vartheta = \{P_{AaAb}, i_{P_{AaAb}}, e_{AaAb}, a_{AaAb}, \omega_{AaAb}, \Omega_{AaAb}, i_{AaAb}, q_{AaAb}\} \cup \{P_{AB}, t_{P_{AB}}, e_{AB}, \omega_{AB}, \Omega_{AB}, i_{AB}, q_{AB}, \pi, V_0\}^6$. We use uniform priors on all the orbital parameters according to its physical valid range (e.g., $q_{AaAb} \sim U_{[0,1]}$), except on the system parallax, where we use a Gaussian prior $\pi \sim N(29.408, 0.027^2)$ mas (according to the value and standard deviation reported by Gaia DR3, Gaia Collaboration 2016b). Based on the spectral types reported for the primary and secondary stellar components reported in Table 1, we also used truncated Gaussian priors on the stellar masses $M_{Aa} \sim N(1.245, 0.067^2) M_\odot$ and $M_B \sim N(0.485, 0.059^2) M_\odot$ within the interval $[\mu -$

<table><caption>**Table 3.** Relative astrometry of HD 196885 B relative to A extracted from NaCo (see Chauvin et al. 2011), SOAR, IRDIS, and GRAVITY (referred as GRAV. in the table) observations.</caption><thead><tr><th>Inst.</th><th>MJD</th><th>Δα<br>(mas)</th><th>Δδ<br>(mas)</th><th>σ<sub>Δα</sub><br>(mas)</th><th>σ<sub>Δα</sub><br>(mas)</th><th>ρ</th><th></th></tr></thead><tbody><tr><td>SOAR</td><td>51826.186</td><td>650.0</td><td>135.8</td><td>5.0</td><td>5.0</td><td></td><td></td></tr><tr><td>SOAR</td><td>52459.690</td><td>709.7</td><td>143.1</td><td>5.0</td><td>5.0</td><td></td><td></td></tr><tr><td>NaCo</td><td>53583.238</td><td>658.5</td><td>273.6</td><td>1.7</td><td>1.6</td><td></td><td></td></tr><tr><td>NaCo</td><td>53974.086</td><td>649.8</td><td>292.6</td><td>1.9</td><td>2.4</td><td></td><td></td></tr><tr><td>NaCo</td><td>54337.053</td><td>640.2</td><td>309.7</td><td>3.1</td><td>3.2</td><td></td><td></td></tr><tr><td>NaCo</td><td>54645.259</td><td>630.9</td><td>326.3</td><td>3.1</td><td>3.2</td><td></td><td></td></tr><tr><td>NaCo</td><td>55070.108</td><td>614.5</td><td>342.1</td><td>3.0</td><td>2.9</td><td></td><td></td></tr><tr><td>IRDIS</td><td>57178.378</td><td>497.0</td><td>405.0</td><td>2.0</td><td>2.0</td><td></td><td></td></tr><tr><td>SOAR</td><td>57291.493</td><td>491.0</td><td>407.7</td><td>3.0</td><td>3.0</td><td></td><td></td></tr><tr><td>IRDIS</td><td>58026.077</td><td>433.0</td><td>422.0</td><td>2.0</td><td>2.0</td><td></td><td></td></tr><tr><td>SOAR</td><td>58265.219</td><td>417.8</td><td>420.7</td><td>3.0</td><td>3.0</td><td></td><td></td></tr><tr><td>SOAR</td><td>58680.881</td><td>383.1</td><td>420.1</td><td>3.0</td><td>3.0</td><td></td><td></td></tr><tr><td>GRAV.</td><td>58710.188</td><td>377.115</td><td>425.575</td><td>0.017</td><td>0.037</td><td>-0.61</td><td></td></tr><tr><td>GRAV.</td><td>59363.331</td><td>320.243</td><td>428.407</td><td>0.015</td><td>0.043</td><td>-0.74</td><td></td></tr><tr><td>GRAV.</td><td>59450.195</td><td>312.573</td><td>428.333</td><td>0.034</td><td>0.049</td><td>-0.96</td><td></td></tr><tr><td>SOAR</td><td>59489.151</td><td>385.2</td><td>3.0</td><td>3.0</td><td>3.0</td><td></td><td></td></tr><tr><td>GRAV.</td><td>59507.082</td><td>307.471</td><td>428.540</td><td>0.038</td><td>0.037</td><td>-0.65</td><td></td></tr></tbody></table>

**Notes.** For GRAVITY, the Pearson's coefficient ρ quantify the correlation between the Δα and Δδ uncertainties.

$2\sigma, \mu + 2\sigma$ (since it accounts ~ 95% of the Gaussian probability density function), using the recent calibration provided by Abushattal et al. (2020)<fnref n="7" />. We run 10000 iterations on each of four independent Markov chains, discarding the first half for the warm-up phase. Each chain was initialized from an starting point determined through the quasi-Newton method L-BFGS (Liu & Nocedal 1989). The code was implemented using the probabilistic programming language Stan (Carpenter et al. 2017).

<footnote n="6">⁶ The semi-major axis of the outer orbit $a_{AB}$ is not part of the parameter space to be sampled since it is deterministically computed using the equivalence between the total mass of the inner system AaAb and the mass of the primary object of the outer system AB, i.e., $M_{Aa} + M_{Ab} = M_A$.</footnote>

<footnote n="7">⁷ Note that the masses $M_{Aa}$, $M_{Ab}$ and $M_B$ are not in the set of parameters to sample since they are generated from $\vartheta$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** *Top-left:* astrometric position of HD 196885 B relative to A combining SOAR, NaCo, SPHERE and GRAVITY observations from August 2005 to August, 2021. Specific zoomed plots on the NaCo astrometric measurements obtained between August, 2005 and August, 2009 (middle-center), the SPHERE and SOAR ones from June, 2015 to May, 2018 (top-right), and the GRAVITY ones from May, 2019 to August, 2021 (middle-right) are shown. *Bottom-left:* radial velocity measurements of HD 196885 A from CORAVEL, ELODIE, CORALIE, Lick and SOPHIE covering 40 yr of monitoring. *Top-right:* zoom on the radial velocity signal induced by HD 196885 Ab. The results of the orbital fitting exploration are reported with the thick blue curves.</figcaption></figure>

## 4.2. A massive Jupiter on a mutually inclined orbit relative to the stellar companion?

The results of the orbital fitting exploration led to the identification of nine for families of solutions (considering a truncated Gaussian prior for the mass of HD 196885 A) that are reported in Table A.1. Their log-posterior density (up to a constant) distributions is also shown in Fig. A.1. In Table 4 and in Fig. 3 (in blue), we report the most likely solutions of the family Sol. 2 which are compared with the observed relative astrometric and radial

velocity measurements of the system. The posterior solutions of each orbital and physical parameters of the exoplanetary companion HD 196885 Ab, and the ones of the stellar companion B are shown in Fig. 4.

The long-term monitoring over almost 40 yr combined with 15 yr of direct imaging set a relatively tight constraints on the orbital properties of the stellar companion. The solutions are very similar to the ones already derived by Chauvin et al. (2011) indicating a moderately eccentric ($e_{AB}$ = 0.417) orbit, a
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Marginal posterior distributions of the physical and orbital parameters of HD 196885 Ab including the period, periastron passage, eccentricity, the semi-major axis... (two first rows). On the two last rows are reported the marginal posterior distributions for HD 196885 B. The maximum a posteriori (the most probable sample of the posterior distribution) are reported in red.</figcaption></figure>

<table><caption>**Table 4.** Summary of the orbital parameters derived for HD 196885 Ab and HD 196885 B: the semi-major axis (a), the eccentricity (e), the inclination (i), the longitude of ascending node (Ω), the argument of periastron (ω), the time of periastron passage (t<sub>p</sub>).</caption><thead><tr><th>Parameters</th><th>Ab</th><th>B</th></tr></thead><tbody><tr><td>a (au)</td><td>2.383<sup>2.385</sup><sub>2.379</sub></td><td>19.778<sup>19.886</sup><sub>19.759</sub></td></tr><tr><td>P (yr)</td><td>3.485<sup>3.486</sup><sub>3.469</sub></td><td>69.045<sup>69.578</sup><sub>68.934</sub></td></tr><tr><td>e</td><td>0.444<sup>0.457</sup><sub>0.439</sub></td><td>0.417<sup>0.418</sup><sub>0.413</sub></td></tr><tr><td>i (°)</td><td>143.041<sup>149.613</sup><sub>138.459</sub></td><td>120.427<sup>120.407</sup><sub>120.183</sub></td></tr><tr><td>Ω (°)</td><td>91.570<sup>116.435</sup><sub>103.696</sub></td><td>79.150<sup>79.138</sup><sub>78.974</sub></td></tr><tr><td>ω (°)</td><td>150.166<sup>159.668</sup><sub>153.034</sub></td><td>231.464<sup>232.036</sup><sub>231.503</sub></td></tr><tr><td>t<sub>p</sub> (yr AD)</td><td>1982.517<sup>1982.732</sup><sub>1982.525</sub></td><td>1982.886<sup>1982.985</sup><sub>1982.832</sub></td></tr><tr><td>M<sub>Aa</sub> (M<sub>☉</sub>)</td><td>1.111<sup>1.114</sup><sub>1.112</sub></td><td></td></tr><tr><td>M<sub>Ab</sub> (M<sub>Jup</sub>)</td><td>3.394<sup>4.096</sup><sub>3.130</sub></td><td></td></tr><tr><td>M<sub>B</sub> (M<sub>☉</sub>)</td><td></td><td>0.508<sup>0.508</sup><sub>0.506</sub></td></tr><tr><td>Φ (°)</td><td></td><td>24.358<sup>37.747</sup><sub>27.596</sub></td></tr></tbody></table>

**Notes.** We report the maximum a posteriori estimates (the most probable solution), with the first and third quartiles (as an uncertainty measure) as subscripts and superscripts, respectively. Note that the MAP values are not necessarily equal to the second quartile (the median).

70 yr period, an inclination of 120.43 deg, and a time to periastron in 1982.87 yr. However, the addition of SOAR, SPHERE, and GRAVITY astrometric measurements from the last decade enables us to nail down previous uncertainties.

For the planet HD 196885 Ab, the gain of almost two orders of magnitude in astrometric accuracy obtained by GRAVITY at VLTI (30–50 µas precision) compared to what was regularly achieved by NaCo and SPHERE at VLT (1–2 mas precision) completely opens a new window to detect and characterize the periodic wobble caused by the S-type planet on the angular separation between both stellar components A and B. The four GRAVITY measurements obtained in August 2019, May 2021, August 2021, and October 2021 covering less than 60% of the orbital period of HD 196885 Ab, allow to unambiguously reject all families of solutions associated with stellar and brown dwarf masses for Ab, and to pinpoint one specific family of retrograde solutions of Jupiter-like planets with low masses ($M_{Ab} = 3.39 M_{Jup}$), with eccentricity ($e_{AaAb} = 0.44$), and inclination close to 143.04 deg with a moderate mutual inclination with B ($Φ = 24.36$ deg).

We note that other islands of solutions remain possible in the planetary mass regime but for different orbital configurations: retrograde and prograde solutions with higher mutual inclination ($Φ = 60–80$ deg or $Φ = 110–120$ deg, respectively) for low masses for Ab ($M_{Ab} = 2–4 M_{Jup}$) and moderate eccentricities, as well as one family of prograde solution for higher masses ($M_{Ab} = 12–14 M_{Jup}$), moderate eccentricity and mutual inclination ($Φ = 130$ deg).
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Dynamical evolution of HD 196885 system as computed with the HJS integrator starting from the orbital solution of Table 4 over $10^8$ yr. *Left-right*: semi-major axes, eccentricities, and mutual inclination. In the *left* and *central* plots, the orange curves stand for the inner orbit and the blue one for the outer orbit, respectively.</figcaption></figure>

**5. Origin and fate**

Our orbital determination clearly suggests that the 3-body system HD 196885 is presumably non-coplanar. Figure 4 shows the current favored family of retrograde solutions for which mutual inclination $\Phi$ lies between typical values of 27.6 deg (first quartile) and 37.8 deg (third quartile; note that the MAP value is equal to 24.4 deg). This non negligible inclination makes it a good candidate for a von Zipele Kozai-Lidov mechanism (or resonance, Kozai 1962) mainly affecting the inner orbit. This dynamical behavior is common in non-planar multiple systems (Beust & Dutrey 2006; Krymolowski & Mazeh 1999; Ford et al. 2000). It is characterized by coupled large amplitude eccentricity and inclination oscillations mainly affecting the inner orbit. Regarding the relatively large individual masses involved here, this could lead the system to instability, which justifies a dedicated stability study.

To do it, we used the HJS integrator (Beust 2003), a symplectic package specifically designed to study hierarchical stellar systems. The same package was successfully indeed used by Beust & Dutrey (2006) to investigate the multiple system GG Tauri. Here we performed an integration over $10^8$ yr starting from the orbital solution from Table 4. The time step was fixed to 0.05 yr, which is a very conservative assumption, as this represents 1/70 of the smaller orbital period (AaAb orbit), while tests by Beust (2003) show that with only 1/20, the total energy is preserved to $10^{-6}$ relative fraction. But as we fear that the system might be close to instability given the high individual masses of the various bodies involved, this appeared to be a reasonable assumption.

Figure 5 shows the result of this computation, taking the orbital fitting solution from Table 4 as starting point. The evolution is shown over $10^8$ yr. This configuration appears stable over the whole integration, as can be seen from the stability of the semi-major axes. Nonetheless the eccentricity and mutual inclination plots suggest von Zipele Kozai Lidov cycles that affect the inner orbit, while the outer orbit is only little affected. This von Zipele Kozai Lidov resonance remains moderate here, due to the rather small mutual inclination. Hence the stability of the triple system is not affected.

The formation mechanism of such a system is unclear. We currently do not know if planets in such tight binaries follow similar or different formation processes than the ones orbiting single stars. For instance, Cadman et al. (2022) have investigated the role of binary companions triggering fragmentation in self-gravitating discs. They find that binarity triggers disk instability for intermediate-close (100–400 au) binaries while the presence of a binary companion within 100 au has a strong negative role against the occurrence of instability. For core accretion, in the

specific case of HD 196885 b, Thebault (2011) conclude that the perturbed dynamical environment is strongly hostile to planetesimal accretion in the region where the planet is located today (note that this study considered a coplanar case). Finally, Marzari & Barbieri (2007) propose that the present dynamical configuration may not reflect the formation process since the binary orbit may have changed in the past once the planet formed. Disrupted triple system or single stellar encounter could be alternative processes. Consequently, this formation and dynamical picture must be further re-investigated in the light of upcoming astrometric observations and more detailed simulations. First, the solution depicted in Table 4 may not be unique, and a dedicated follow-up with GRAVITY should enable to rapidly secure the most probable family of solutions. Due to the complexity of the parameter space, other configurations might fit the data as well. Second, other effects could need to be taken into account such as a combination of von Zipele Kozai Lidov cycles and tides (e.g. Fabrycky & Tremaine 2007; Correia et al. 2011; Beust et al. 2012).

**6. Conclusion**

In this study, we have revisited the global architecture of the extreme planetary system HD 196885. This system is composed of a tight F8V-M1V binary hosting a S-type circumprimary planet around the A component. Based on the combination of new measurements in radial velocity, speckle interferometry, high-contrast imaging and high-precision astrometry with interferometry, we could narrow down the orbital properties of both HD 196885 B and HD 196885 Ab. The binary companion orbits on an inclined and retrograde plan ($i_{AB} = 120.43$ deg) with a semi-major axis of 19.78 au, and an eccentricity of 0.415 tightly constrained by almost 40 yr and 20 yr of radial and astrometric monitoring, respectively. For the radial velocity planet HD 196885 Ab, we unambiguously confirm its mass ($M_{Ab} = 3.39 M_{Jup}$) and planetary nature, rejecting all families of pole-on solutions corresponding to stellar and brown dwarf masses. The most favored island of orbital solutions found is associated with moderate eccentricity ($e_{AaAb} = 0.46$), and moderate-inclination close to 143.04 deg. This results points toward a moderate mutual inclination ($\Phi = 24.36$ deg) between the orbital plans of B and Ab. This configuration seems to be dynamically stable over time, although the eccentricity and mutual inclination variations clearly reveal the presence of moderate von Zipele Kozai Lidov cycles that may affect the inner planet. Further observations in the coming years will allow to reject remaining islands of orbital solutions for Ab. Further dynamical simulation should also enable to incorporate additional effects as tides to assess the inner orbit dynamical evolution and stability.
**Acknowledgements.** We would like to thank the staff of ESO VLT/I for their support at the telescope at Paranal and La Silla, and the preparation of the observation at Garching. This publication made use of the SIMBAD and VizieR database operated at the CDS, Strasbourg, France. 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. This work is partly based on data products produced at the SPHERE Data Centre hosted at OSUG/IPAG, Grenoble. We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process. RAM and MV acknowledge partial support from FONDECYT/ANID Grant No. 1190038 and to the Chilean National Time Allocation Committee for their continuous support of our SOAR Speckle survey of southern binaries and multiple systems through programs CN2018A-1, CN2019A-2, CN2019B-13, and CN2020A-19, CN2020B-10, and CN2021B-17. CV and MV acknowledge support from the visiting program of the CNRS French-Chilean Laboratory for Astrophysics (IRL-3386). VF acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grant No. 80NSSC21K0394 (PI: S. Ertel). Based in part on radial velocities obtained with the SOPHIE spectrograph mounted on the 1.93 m telescope at Observatoire de Haute-Provence. We warmly thank the OHP staff for their support on the observations. I.B. received funding from the French Programme National de Planétologie (PNP) of CNRS (INSU). A.C. acknowledges support from CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), PHOBOS (POCI-01-0145-FEDER-029932), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal. JSJ greatly acknowledges support by FONDECYT grant 1201371 and from the ANID BASAL projects ACE210002 and FB210003. A.G. acknowledges support from the ANID-ALMA fund no. ASTRO20-0059.

## References

Abushattal, A. A., Docobo, J. A., & Campo, P. P. 2020, AJ, 159, 28

Baranec, C., Ziegler, C., Law, N. M., et al. 2016, AJ, 152, 18

Bazsó, Á., Pilat-Lohinger, E., Eggl, S., et al. 2017, MNRAS, 466, 1555

Benedict, G. F., Harrison, T. E., Endl, M., & Torres, G. 2018, RNAAS, 2, 7

Beust, H. 2003, A&A, 400, 1129

Beust, H., & Dutrey, A. 2006, A&A, 446, 137

Beust, H., Bonfils, X., Montagnier, G., Delfosse, X., & Forveille, T. 2012, A&A, 545, A88

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Bonavita, M., & Desidera, S. 2007, A&A, 468, 721

Bonavita, M., & Desidera, S. 2020, Galaxies, 8, 16

Bouchy, F., Hébrard, G., Udry, S., et al. 2009, A&A, 505, 853

Bouchy, F., Díaz, R. F., Hébrard, G., et al. 2013, A&A, 549, A49

Cadman, J., Hall, C., Fontanive, C., & Rice, K. 2022, MNRAS, 511, 457

Carpenter, B., Gelman, A., Hoffman, M. D., et al. 2017, J. Stat. Softw., 76, 1

Chauvin, G., Lagrange, A. M., Udry, S., et al. 2006, A&A, 456, 1165

Chauvin, G., Lagrange, A. M., Udry, S., & Mayor, M. 2007, A&A, 475, 723

Chauvin, G., Beust, H., Lagrange, A. M., & Eggenberger, A. 2011, A&A, 528, A8

Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 331

Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, SPIE Conf. Ser., 7014

Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271

Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105

Desidera, S., & Barbieri, M. 2007, A&A, 462, 345

Desidera, S., Gratton, R., Endl, M., et al. 2010, Astrophys. Space Sci. Lib. 366, 105

Desidera, S., Chauvin, G., Bonavita, M., et al. 2021, A&A, 651, A70

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Proc. SPIE, 7014, 70143L

Duchêne, G. 2010, ApJ, 709, L114

Dvorak, R. 1982, Oesterreich. Akad. Wiss. Math. Naturw. Klasse Sitzungsberichte Abteilung, 191, 423

Eggenberger, A., & Udry, S. 2007, arXiv e-prints [arXiv:0705.3173]

Eggenberger, A., Udry, S., Mayor, M., et al. 2004, ASP Conf. Ser., 321, 93

Eggenberger, A., Udry, S., Chauvin, G., et al. 2007, A&A, 474, 273

Eggenberger, A., Udry, S., Chauvin, G., et al. 2008, ASP Conf. Ser., 398, 179

Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298

Fischer, D., Driscoll, P., Isaacson, H., et al. 2009, ApJ, 703, 1545

Fontanive, C., & Bardalez Gagliuffi, D. 2021, Front. Astron. Space Sci., 8, 16

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

Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385

Fragner, M. M., Nelson, R. P., & Kley, W. 2011, A&A, 528, A40

Gaia Collaboration (Brown, A. G. A., et al.) 2016a, A&A, 595, A2

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

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Ginski, C., Mugrauer, M., Adam, C., Vogt, N., & van Holstein, R. G. 2021, A&A, 649, A156

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

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

Gregory, P. 2005, ApJ, 631, 1198

Gregory, P. C. 2011, MNRAS, 410, 94

Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383

Heidari, N., Boisse, I., Orell-Miquel, J., et al. 2022, A&A, 658, A176

Hoffman, M. D., Gelman, A., et al. 2014, J. Mach. Learn. Res., 15, 1593

Hou, F., Goodman, J., Hogg, D. W., Weare, J., & Schwab, C. 2012, ApJ, 745, 198

Jang-Condell, H. 2015, ApJ, 799, 147

Jenkins, J. S., Murgas, F., Rojo, P., et al. 2011, A&A, 531, A8

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

Konacki, M. 2005, AAS Meet. Abstr., 207, 84.02

Kozai, Y. 1962, AJ, 67, 591

Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720

Langlois, M., Gratton, R., Lagrange, A.-M., et al. 2021, A&A, 651, A71

Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, SPIE Conf. Ser., 9146, 91462D

Law, N. M., Morton, T., Baranec, C., et al. 2014, ApJ, 791, 35

Liu, D. C., & Nocedal, J. 1989, Math. Program., 45, 503

Luhman, K. L., & Jayawardhana, R. 2002, ApJ, 566, 1132

Marzari, F., & Barbieri, M. 2007, A&A, 467, 347

Michel, K. U., & Mugrauer, M. 2021, Front. Astron. Space Sci., 8, 14

Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2008, ASP Conf. Ser., 398, 235

Mugrauer, M. 2019, MNRAS, 490, 5088

Mugrauer, M., & Ginski, C. 2015, MNRAS, 450, 3127

Nelson, B., Ford, E. B., & Payne, M. J. 2013, ApJS, 210, 11

Paardekooper, S. J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973

Patience, J., White, R. J., Ghez, A. M., et al. 2002, ApJ, 581, 654

Queloz, D., Mayor, M., Weber, L., et al. 2000, A&A, 354, 99

Raghavan, D., Henry, T. J., Mason, B. D., et al. 2006, ApJ, 646, 523

Silsbee, K., & Rafikov, R. R. 2021, A&A, 652, A104

Tamuz, O., Ségransan, D., Udry, S., et al. 2008, A&A, 480, L33

Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29

Thebault, P., & Haghighipour, N. 2015, in Planetary Exploration and Science: Recent Results and Advances, 309

Thebault, P., Marzari, F., & Scholl, H. 2006, Icarus, 183, 193

Tokovinin, A. 2018, PASP, 130, 035002

Tokovinin, A., & Cantarutti, R. 2008, PASP, 120, 170

Tokovinin, A., Mason, B. D., & Hartkopf, W. I. 2010, AJ, 139, 743

Tokovinin, A., Mason, B. D., Mendez, R. A., et al. 2021, AJ, 162, 41

Toyota, E., Itoh, Y., Ishiguma, S., et al. 2009, PASJ, 61, 19

Videla, M., Mendez, R. A., Clavería, R. M., Silva, J. F., & Orchard, M. E. 2022, AJ, 163, 220

Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, 407, 71

Villegas, C., Mendez, R. A., Silva, J. F., & Orchard, M. E. 2021, PASP, 133, 074501

Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2004, ApJS, 152, 261

Ziegler, C., Tokovinin, A., Latiolais, M., et al. 2021, AJ, 162, 192

Zucker, S., & Mazeh, T. 2002, ApJ, 568, L113

<footnote n="1">¹ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386 and Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="2">² Laboratoire J.-L. Lagrange, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, 06304 Nice, France e-mail: gael.chauvin@oca.eu</footnote>

<footnote n="3">³ Department of Electrical Engineering, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Tupper 2007, Santiago, Chile</footnote>

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

<footnote n="5">⁵ Department of Astronomy, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="6">⁶ CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal</footnote>

<footnote n="7">⁷ IMCCE, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>
⁸ LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France

⁹ Cerro Tololo Inter-American Observatory, Casilla 603, La Serena, Chile

¹⁰ Geneva Observatory, University of Geneva, 51 ch. Pegasi, 1290 Versoix, Switzerland

¹¹ Aix-Marseille Univ, CNRS, CNES, LAM, 13000 Marseille, France

¹² School of Physical Sciences, Faculty of Science, Technology, Engineering and Mathematics, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK

¹³ INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

¹⁴ Steward Observatory, Department of Astronomy, University of Arizona, 933 N. Cherry Ave, Tucson, AZ 85721, USA

¹⁵ Universidad de Concepción, Departamento de Astronomía, Casilla 160-C, Concepción, Chile

¹⁶ European Southern Observatory, Alonso de Cordova 3107, Casilla 19001, Vitacura, Santiago, Chile

¹⁷ Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

¹⁸ Centro de Astrofísica y Tecnologías A fines (CATA), Casilla 36-D, Santiago, Chile
## Appendix A: Families of three-bodies fitting solutions using only astrometric datasets from SOAR, NaCo, SPHERE and GRAVITY

In Fig. A.1 are reported the log posterior values of the different orbital fitting solutions found for the truncated gaussian prior for the mass of the primary. In Table A.1 are summarized the orbital solutions found for these different families.

<figure><img src="image_7.png" /><figcaption>Fig. A.1. Log-posterior density (up to a constant) distributions of 9 computed inferences considering truncated Gaussian priors on the stellar masses.</figcaption></figure>

<table><caption>Table A.1. Families of solutions for the truncated gaussian prior for the mass of HD 196885 A.</caption><thead><tr><th>Sol. #</th><th>M<sub>Aa</sub> (M<sub>☉</sub>)</th><th>M<sub>Ab</sub> (M<sub>Jup</sub>)</th><th>e<sub>AaAb</sub></th><th>i<sub>AaAb</sub> (°)</th><th>M<sub>B</sub> (M<sub>☉</sub>)</th><th>e<sub>AB</sub></th><th>i<sub>AB</sub> (°)</th><th>Φ (°)</th></tr></thead><tbody><tr><td>1</td><td>1.113<sup>1.113</sup><sub>1.111</sub></td><td>3.357<sup>8.793</sup><sub>3.375</sub></td><td>0.038<sup>0.117</sup><sub>0.060</sub></td><td>31.107<sup>30.922</sup><sub>10.883</sub></td><td>0.526<sup>0.530</sup><sub>0.528</sub></td><td>0.406<sup>0.400</sup><sub>0.385</sub></td><td>119.964<sup>119.727</sup><sub>119.073</sub></td><td>117.138<sup>119.184</sup><sub>115.951</sub></td></tr><tr><td>2</td><td>1.111<sup>1.114</sup><sub>1.112</sub></td><td>3.394<sup>4.096</sup><sub>3.130</sub></td><td>0.444<sup>0.457</sup><sub>0.439</sub></td><td>143.041<sup>149.613</sup><sub>138.459</sub></td><td>0.508<sup>0.508</sup><sub>0.506</sub></td><td>0.417<sup>0.418</sup><sub>0.413</sub></td><td>120.427<sup>120.407</sup><sub>120.183</sub></td><td>24.358<sup>37.747</sup><sub>27.596</sub></td></tr><tr><td>3</td><td>1.112<sup>1.114</sup><sub>1.112</sub></td><td>4.877<sup>4.274</sup><sub>2.501</sub></td><td>0.430<sup>0.443</sup><sub>0.428</sub></td><td>155.006<sup>151.312</sup><sub>125.213</sub></td><td>0.511<sup>0.511</sup><sub>0.509</sub></td><td>0.417<sup>0.419</sup><sub>0.413</sub></td><td>120.680<sup>120.669</sup><sub>120.422</sub></td><td>58.363<sup>66.608</sup><sub>58.865</sub></td></tr><tr><td>4</td><td>1.111<sup>1.112</sup><sub>1.111</sub></td><td>4.647<sup>9.764</sup><sub>5.481</sub></td><td>0.149<sup>0.017</sup><sub>0.011</sub></td><td>152.503<sup>168.088</sup><sub>158.225</sub></td><td>0.523<sup>0.534</sup><sub>0.525</sub></td><td>0.410<sup>0.418</sup><sub>0.410</sub></td><td>120.531<sup>121.003</sup><sub>120.489</sub></td><td>65.833<sup>69.971</sup><sub>65.140</sub></td></tr><tr><td>5</td><td>1.111<sup>1.111</sup><sub>1.111</sub></td><td>7.144<sup>7.202</sup><sub>7.125</sub></td><td>0.402<sup>0.403</sup><sub>0.401</sub></td><td>14.256<sup>14.378</sup><sub>14.266</sub></td><td>0.537<sup>0.536</sup><sub>0.534</sub></td><td>0.382<sup>0.387</sup><sub>0.384</sub></td><td>119.691<sup>119.800</sup><sub>119.658</sub></td><td>131.891<sup>132.162</sup><sub>131.863</sub></td></tr><tr><td>6</td><td>1.111<sup>1.112</sup><sub>1.111</sub></td><td>12.129<sup>12.744</sup><sub>12.349</sub></td><td>0.319<sup>0.342</sup><sub>0.326</sub></td><td>8.872<sup>8.770</sup><sub>8.515</sub></td><td>0.550<sup>0.552</sup><sub>0.550</sub></td><td>0.369<sup>0.371</sup><sub>0.369</sub></td><td>119.032<sup>119.186</sup><sub>119.010</sub></td><td>125.774<sup>126.198</sup><sub>125.910</sub></td></tr><tr><td>7</td><td>1.111<sup>1.112</sup><sub>1.111</sub></td><td>3.497<sup>3.639</sup><sub>2.372</sub></td><td>0.453<sup>0.481</sup><sub>0.451</sub></td><td>144.252<sup>145.394</sup><sub>55.234</sub></td><td>0.538<sup>0.542</sup><sub>0.538</sub></td><td>0.385<sup>0.389</sup><sub>0.382</sub></td><td>119.218<sup>119.470</sup><sub>119.145</sub></td><td>25.696<sup>96.359</sup><sub>54.869</sub></td></tr><tr><td>8</td><td>1.111<sup>1.112</sup><sub>1.111</sub></td><td>8.922<sup>10.630</sup><sub>9.901</sub></td><td>0.459<sup>0.477</sup><sub>0.458</sub></td><td>13.124<sup>11.745</sup><sub>10.861</sub></td><td>0.540<sup>0.541</sup><sub>0.539</sub></td><td>0.401<sup>0.414</sup><sub>0.404</sub></td><td>119.544<sup>120.042</sup><sub>119.649</sub></td><td>111.104<sup>114.008</sup><sub>112.750</sub></td></tr><tr><td>9</td><td>1.111<sup>1.112</sup><sub>1.111</sub></td><td>3.657<sup>4.017</sup><sub>3.654</sub></td><td>0.549<sup>0.569</sup><sub>0.439</sub></td><td>151.370<sup>155.523</sup><sub>151.634</sub></td><td>0.539<sup>0.543</sup><sub>0.539</sub></td><td>0.412<sup>0.412</sup><sub>0.398</sub></td><td>120.230<sup>120.242</sup><sub>119.703</sub></td><td>88.322<sup>88.129</sup><sub>84.440</sub></td></tr></tbody></table>