# An imaged 15 M Jup companion within a hierarchical quadruple system

A. Chomez, V. Squicciarini, A.-M. Lagrange, P. Delorme, G. Viswanath, M.

Janson, O. Flasseur, G. Chauvin, M. Langlois, P. Rubini, et al.

► To cite this version:

A. Chomez, V. Squicciarini, A.-M. Lagrange, P. Delorme, G. Viswanath, et al. An imaged 15 M Jup companion within a hierarchical quadruple system. Astronomy & Astrophysics - A&A, 2023, 676, pp.L10. 10.1051/0004-6361/202347044. hal-04769778

HAL Id: hal-04769778

https://hal.science/hal-04769778v1

Submitted on 6 Nov 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

# An imaged 15 $M_{Jup}$ companion within a hierarchical quadruple system*

A. Chomez<sup>1,2</sup>, V. Squicciarini<sup>1,3</sup>, A.-M. Lagrange<sup>1,2</sup>, P. Delorme<sup>2</sup>, G. Viswanath<sup>4</sup>, M. Janson<sup>4</sup>, O. Flasseur<sup>5</sup>, G. Chauvin<sup>6</sup>, M. Langlois<sup>5</sup>, P. Rubini<sup>7</sup>, S. Bergeon<sup>2</sup>, D. Albert<sup>16</sup>, M. Bonnefoy<sup>2</sup>, S. Desidera<sup>3</sup>, N. Engler<sup>8</sup>, R. Gratton<sup>3</sup>, T. Henning<sup>9</sup>, E. E. Mamajek<sup>10</sup>, G.-D. Marleau<sup>11,12,13,9</sup>, M. R. Meyer<sup>14</sup>, S. Reffert<sup>15</sup>, S. C. Ringqvist<sup>4</sup>, and M. Samland<sup>4</sup>

*(Affiliations can be found after the references)*

Received 29 May 2023 / Accepted 30 June 2023

## ABSTRACT

**Context.** Since 2019, the direct imaging B-star Exoplanet Abundance STudy (BEAST) at SPHERE@VLT has been scanning the surroundings of young B-type stars in order to ascertain the ultimate frontiers of giant planet formation. Recently, the $17_{-4}^{+3}$ Myr HIP 81208 was found to host a close-in (~50 au) brown dwarf and a wider (~230 au) late M star around the central $2.6 M_{\odot}$ primary.

**Aims.** Alongside the continuation of the survey, we are undertaking a complete reanalysis of archival data aimed at improving detection performances so as to uncover additional low-mass companions.

**Methods.** We present here a new reduction of the observations of HIP 81208 using the patch covariance algorithm (PACO), a recent and powerful algorithm dedicated to processing high-contrast imaging datasets, as well as more classical algorithms and a dedicated point spread function subtraction approach. The combination of different techniques allowed for a reliable extraction of astrometric and photometric parameters.

**Results.** A previously undetected source was recovered at a short separation from the C component of the system. Proper motion analysis provided robust evidence for the gravitational bond of the object to HIP 81208 C. Orbiting C at a distance of ~20 au, this 15 $M_{Jup}$ brown dwarf becomes the fourth object of the hierarchical HIP 81208 system.

**Conclusions.** Among the several BEAST stars which are being found to host substellar companions, HIP 81208 stands out as a particularly striking system. As the first stellar binary system with substellar companions around each component ever found by direct imaging, it yields exquisite opportunities for thorough formation and dynamical follow-up studies.

**Key words.** techniques: high angular resolution – planetary systems – brown dwarfs – stars: individual: HIP 81208 – planets and satellites: detection

## 1. Introduction

The formation of planets in binary systems, and chiefly the tightest (≤50 au) ones, is a vibrant subject in exoplanetology. Binary systems are indeed complex environments from a dynamical point of view, severely affecting the size of protoplanetary disks and their capability to either form massive enough cores to undergo runaway gas accretion that is core accretion (CA; Bodenheimer & Pollack 1986) or induce low enough Toomre Q values to trigger gravitational instability (GI; Boss 1997). Whether substellar companions can form critically depends on the stars properties, their physical separations, and the disk initial properties (Bouwman et al. 2006; Jang-Condell 2015; Silsbee & Rafikov 2021).

From an observational standpoint, ~300 S-type substellar companions (companions orbiting one component of a binary system) within binary systems are known to date (Fontanive & Bardalez Gagliuffi 2021; Chauvin et al. 2023) – with their frequency being anticorrelated with binary separation (Wang et al. 2014) – as well as triple and higher-order planet-hosting systems in strongly hierarchical configurations (Roberts 2015; Cuntz et al. 2022). Indirect techniques have identified a few systems where both components host substellar companions (see, e.g., Desidera et al. 2014; Lissauer et al. 2014; Udry et al. 2019);

notably, HD 41004 stands out due to its close A-B separation (~23 au): with a $m \sin i = 18 M_{Jup}$ brown dwarf orbiting at 0.017 au from component B ($0.4 M_{\odot}$) and a $m \sin i = 2.5 M_{Jup}$ companion around component A ($0.7 M_{\odot}$) on an orbit with semi-major axis $a = 1.6$ au and a large eccentricity $e = 0.4$ (Zucker et al. 2003, 2004).

In the course of a new analysis of archival data obtained through the Spectro-Polarimetric High contrast imager for Exoplanets REsearch (SPHERE; Beuzit et al. 2019), we detected a new companion in the young ($17_{-4}^{+3}$ Myr) HIP 81208 system. HIP 81208 was observed as part of the B-star Exoplanet Abundance STudy (BEAST) survey dedicated to the search for exoplanets around 85 B-type members of the Scorpius Centaurus (Sco-Cen) association (Janson et al. 2021a). Located in the Upper Centaurus-Lupus (UCL) subgroup of Sco-Cen at a distance of $148.7_{-1.3}^{+1.5}$ pc (Gaia Collaboration 2023), it has been recently identified as a triple system, where the A component is a $2.58 \pm 0.06 M_{\odot}$ B9V star, the B component is a $67_{-7}^{+6} M_{Jup}$ brown dwarf orbiting HIP 81208 A at about 50 au, and the C component is a low-mass star of $0.135_{-0.013}^{+0.010} M_{\odot}$ orbiting HIP 81208 A at about 230 au (Viswanath et al. 2023). The newly found companion (hereafter Cb, following naming conventions for hierarchical systems; Hartkopf & Mason 2004; Deeg & Belmonte 2018) is orbiting the C component, making HIP 81208 the first binary system with substellar companions around both components ever discovered through direct imaging (DI).

<footnote>* Based on data obtained with the ESO/VLT SPHERE instrument under programs 1101.C-0258(A/E).</footnote>
<table><caption>**Table 1.** Observation logs for the two epochs.</caption><thead><tr><th></th><th>First epoch (2019-08-06)</th><th>Second epoch (2022-04-05)</th></tr></thead><tbody><tr><td>IRDIS filter</td><td>DB_K12</td><td>DB_K12</td></tr><tr><td>DIT(s) × Nframe</td><td>64 × 48</td><td>64 × 48</td></tr><tr><td>ΔPA (°)</td><td>59.09</td><td>57.33</td></tr><tr><td>Seeing (′′) <sup>(a)</sup></td><td>0.64</td><td>0.57</td></tr><tr><td>Airmass <sup>(a)</sup></td><td>1.020</td><td>1.020</td></tr><tr><td>τ<sub>0</sub> (ms) <sup>(a)</sup></td><td>8.1</td><td>5.3</td></tr><tr><td>Program ID</td><td>1101.C-0258(A)</td><td>1101.C-0258(E)</td></tr></tbody></table>

**Notes.** DIT = detector integration time per frame, ΔPA = amplitude of the parallactic rotation, τ₀ = coherence time. <sup>(a)</sup>Values extracted from the updated DIMM information and averaged over the sequence.

We present SPHERE data and data reduction in Sect. 2; after confirming the bound nature of the companion candidate, we describe its properties in Sect. 3. A discussion on the peculiar properties of this quadruple system follows in Sect. 4.

## 2. Data analysis

### 2.1. SPHERE data

HIP 81208 was observed twice by SPHERE (Beuzit et al. 2019) on August 6, 2019 and on April 5, 2022. Both observations were conducted using the telescope in pupil-stabilized mode. This allows for the use of angular and spectral differential imaging (ASDI; Marois et al. 2006) post-processing techniques. In each case, an unsaturated, noncoronagraphic image – point spread function (PSF) – of the primary was obtained for flux calibration purposes, as well as a coronagraphic exposure with a waffle pattern applied to the mirror (Cantalloube et al. 2019), for centering purposes, before and after the main coronagraphic exposures. The N-ALC-YJH-S coronagraph was used, allowing the infrared dual-band imager and spectrograph (IRDIS; Dohlen et al. 2008) to observe in the K1 and K2 bands while the integral field spectrograph (IFS; Claudi et al. 2008) observed in the YJH bands. Because the source of interest for this Letter is outside the field of view of IFS, only IRDIS data are considered hereafter. Table 1 summarizes the observing conditions as well as the setup for the two observations, with the same already having been used in Viswanath et al. (2023).

### 2.2. Data reduction

Data reduction was performed on the COBREX Data Center, which is a modified and improved server based on the High Contrast Data Center, (HC-DC, formerly SPHERE Data Center; Delorme et al. 2017). The COBREX Data Center aims to improve the detection capabilities with existing SPHERE images by means of the patch covariance (PACO) algorithm. More specifically, we used PACO ASDI (Flasseur et al. 2018, 2020b,a) as well as the No-ADI routine embedded in the SPHERE speckle calibration software (SPECIAL; Galicher et al. 2018) as post-processing algorithms. The prereduction pipeline (i.e., going from raw data to a calibrated 4D datacube) is identical to the one implemented in HC-DC, performing dark, flat, distortion, and bad pixel corrections.

PACO models the noise using a multi-Gaussian model at a local scale on small patches, allowing for a better estimation of the temporal and spectral correlation of the noise. The full details on the improvements of the prereduction pipeline

as well as the optimization regarding the ASDI mode of PACO, and the obtained performances are described in Chomez et al. (2023). PACO provides a contrast gain between one and two magnitudes at all separations as well as reliable and statistically grounded signal-to-noise ratio (S/N) detection and contrast maps in an unsupervised fashion compared to more classical algorithms such as TLOCI-ADI used in SPECIAL (Galicher et al. 2018). Prompted by the results emerging from our PACO and No-ADI reductions (Sect. 2.3), we additionally developed a custom PSF subtraction routine, building a local PSF model to remove C and enhance detection capabilities in its immediate surroundings (see the detailed description in Appendix A).

In regards to true north, the pixel scale, and the pupil offset, we adopted the long-term IRDIS calibrations by Maire et al. (2021): a pixel scale of 12.258 ± 0.004 mas pixel⁻<fnref n="1" /> (K1 band) and 12.253 ± 0.003 mas/pixel (K2 band), a true north orientation of (−1.77 ± 0.04)°, and a pupil offset of (136 ± 0.03)°.

### 2.3. Detected sources

Figure 1 shows the S/N maps of PACO and residual No-ADI maps for both epochs: a new source (hereafter Cb) was detected in close proximity (~120 mas) to C.

PACO detects Cb with a high S/N (26.2 in 2019 and 16.9 in 2022), placing a high confidence level on the detection. Although Cb is visible by eye in the no-ADI reduction, no reliable measurement could be extracted for the source because of strong contamination by C; SPECIAL does not provide tools to deconvolve sources. Conversely, PACO includes a “cleaning” option designed for the case: it removes the contribution of C while characterizing Cb, enabling extraction of a reliable photometry (Flasseur et al. 2020b). We attribute the previous nondetection of Cb to self-subtraction artifacts introduced by TLOCI – the baseline algorithm used to process IRDIS observations in BEAST’s standard reduction pipeline – near C<fnref n="1" />.

Our PSF subtraction routine, specifically designed to investigate the surroundings of C, enabled us to solidly reveal Cb and its first Airy ring on almost 360 degrees (see Fig. 2) at both epochs. Figure A.1 shows the residuals after removing both C and Cb. The highest residuals in both epochs are barely above the local background noise in each of the 48 individual frames, allowing us to robustly determine the position, contrast, and associated uncertainties by deriving the mean and standard deviation of the resulting 48 independent measurements. Notably, the source is not only visible with both algorithms at both epochs, but also in each raw frame before any post-processing algorithm is applied (Fig. B.1); unlike the nearby first Airy ring of C, its separation from C does not vary with wavelength, while its rotation around C during the ADI sequence is consistent with the expected motion of a physical source (see Appendix B), thus ensuring it is not an artifact.

As an additional check, we also characterized B and C, finding results compatible within 1σ to those presented by Viswanath et al. (2023)<fnref n="2" />. Furthermore, besides redetecting all known background sources with similar astrometric and photometric values to the initial analysis, we imaged an additional faint source (CC14) owing to PACO’s deeper sensitivity.

<footnote n="1">¹ Any ADI-based algorithm with a subtraction step (e.g., which does not fit the planet and the systematics simultaneously) also suffers from this self-subtraction effect.</footnote>

<footnote n="2">² We attribute the larger astrometric uncertainties emerging in our analysis to the fact that the previous analysis did not include primary centering uncertainties, which dominate positional uncertainty here.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** S/N maps produced by PACO (top row) and residual maps produced by SPECIAL no-ADI (bottom row) for both epochs (with 2019 on the left and 2022 on the right). The B and C components were retrieved and the reductions unveil the additional Cb companion. The central star (hidden behind the coronagraphic mask) is indicated by the yellow star.</figcaption></figure>

Astrophotometrical results for B, C, Cb, and CC14 are provided in Appendix C.

**3. A low-mass substellar companion around HIP 81208 C**

**3.1. Companion confirmation**

Figure 3 shows the proper motion of B, C, and Cb as opposed to the already known background sources. As for B and C, the motion of Cb is inconsistent (at $6\sigma$) with the observed motion of field interlopers. We anticipate that the non-null relative motion between C and Cb is consistent with orbital motion, as it is subsequently quantified in Sect. 3.3.

Having confirmed the common proper motion of Cb to the already known HIP 81208 A, B, and C, we carefully investigated an alternative possibility: namely, that A-B and C-Cb constitute two independent binary members of UCL projected at a short separation (Appendix D). We therefore assessed the probability that an UCL member, unseen by Gaia (Gaia Collaboration 2023), could end up as an interloper to any BEAST star. A final false alarm probability, referring to at least one detection across the entire current survey (47 stars), of $1.3 \times 10^{-3}$ was obtained (meaning that we expect, on average, one false positive out of

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Zoom on the C-Cb pair after subtracting the PSF of C (with 2019 on the left and 2022 on the right) for the K1 band. The position of C (subtracted) is highlighted by the white cross. We note that Cb is clearly visible after the removal of C.</figcaption></figure>

~40 000 observed B stars), placing a high confidence level on the membership of A (+B) and C (+Cb) to a single quadruple system.

**3.2. Physical properties**

As in Viswanath et al. (2023), an estimate of the photometric mass of the newly discovered object was obtained using the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Astrometric displacements between the two epochs for all known sources (data from Viswanath et al. 2023) with the addition of Cb (blue diamond) and CC14 (red downward-pointing triangle). The expected displacement at the second epoch is indicated by a white and a black star for a static background source and a comoving source with the target, respectively. The mean astrometric shift of background sources (gray crosses) is shown as a thick gray cross. Inset: Relative motion between Cb and C, with C set as the primary; the error bar associated with Cb is small enough to be hidden behind the diamond.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Color magnitude diagram with the three components (B, C, and Cb) represented by the red dots.</figcaption></figure>

MADYS tool (Squicciarini & Bonavita 2022): we combined PACO (K1, K2) contrasts, the 2MASS $K_s$ magnitude of the primary (Cutri et al. 2003), plus the system's color excess ($E(B-V) = 0.011 \pm 0.021$ mag) and age ($17_{-4}^{+3}$ Myr). For the purpose of accounting for theoretical uncertainties on the final estimates, the computation was performed by comparison with two different models suited for the age and mass range of interest: the Ames-Dusty models (Allard et al. 2001) and the BT-Settl models (Allard 2016). The resulting values were averaged to yield a final mass estimate:

$$M_{Cb} = 14.8 \pm 0.4 M_{Jup}. \quad (1)$$

Additional details on the derivation of this estimate and its associated uncertainty are provided in Appendix E. Based on the models, this best-fit mass would correspond to expected 2MASS $H = 10.28 \pm 0.07$ mag and $K_s = 9.74 \pm 0.04$ mag.

In the same fashion, the average effective temperature, surface gravity, and bolometric luminosity returned by this comparison are $T_{eff} = 2050_{-20}^{+35}$ K, $\log g = 4.087_{-0.022}^{+0.011}$, and $\log L/L_{\odot} = -3.31 \pm 0.03$, respectively.

However, we acknowledge that we currently have only two photometric points probing similar wavelengths and that, even in the case of young substellar objects with much more comprehensive data available, systematic errors intrinsic to atmosphere and evolution models might be up to an order of magnitude larger than formal uncertainties (see, e.g., Petrus et al. 2021). These differences can arise, for instance, from uncertainties on the initial entropy after accretion, possible age differences between a planet and host star, as well theoretical difficulties in handling atmospheric dust (see, e.g., Lueber et al. 2023).

Figure 4 displays the position of B, C, and Cb on a color-magnitude diagram, with all matching the M sequence, while

Table 2 reports the main outputs of the astrophotometric characterization of the object.

### 3.3. Orbital properties

We ran the emcee code (Foreman-Mackey et al. 2013) to derive information on the orbital parameters of the companion starting from the astrometry and their best-fit masses of Cb and C. We used their relative astrometry, as measured by PSF subtraction, because it is affected by much smaller uncertainties than absolute astrometry (Appendix A). The sampling tool is based on the emcee 3.0 library, using a mix of custom move functions to alleviate potential multimodality problems and the cyclicity of angular variables. We note that 100 walkers, between 100 000 and 400 000 iterations, and ten temperatures were used. The priors include a uniform log prior for $sma$ ($a \in [0, 80]$ au). The upper value in sma corresponds roughly to 0.3 times the projected separation between A and C, which, following Musielak et al. (2005), should allow for the stability of Cb in the binary system, given the masses of A and C, and assuming a null eccentricity. We nonetheless considered the eccentricities' $e \in [0, 0.4]$ range for the priors<fnref n="3" />.

Given the limited information available, the orbital parameters are poorly constrained (see Fig. F.1 and in Appendix F). The $a$ distribution peaks at 17 au ($T \sim 190$ yr), with a tail extending to more than 40 au ($T \sim 600$ yr). The eccentricity is not constrained. The inclination of its orbital plane is $i = 73 \pm 20^{\circ}$. Figure 5 shows the 1000 best draws from our posterior distributions for Cb determined based on the loglikelihood. We also ran the MCMC on B and C companions, finding $a$, $e$, and $i$ compatible with those found by Viswanath et al. (2023), albeit with larger uncertainties, due to the larger error bars found in the present astrometric measurements.

<footnote n="3">³ We note that assuming a larger range (up to 0.9) does not significantly change the results.</footnote>
<table><caption>**Table 2.** Final astrometry for Cb with respect to A (sep, PA) and C (ΔRA, ΔDec), measured IRDIS contrasts and best-fit values for mass and $T_{\text{eff}}$.</caption><thead><tr><th>Epoch</th><th>Sep (mas)</th><th>PA (°)</th><th>ΔK<sub>1</sub> (mag)</th><th>ΔK<sub>2</sub> (mag)</th><th>ΔRA (mas)</th><th>ΔDec (mas)</th><th>Mass (M<sub>Jup</sub>)</th><th>q<sub>Cb</sub></th><th>T<sub>eff</sub> (K)</th></tr></thead><tbody><tr><td>2019</td><td>1573.4 ± 3.2</td><td>119.9 ± 0.12</td><td>8.93 ± 0.08</td><td>8.59 ± 0.08</td><td>26.4 ± 1.13</td><td>-116.9 ± 0.99</td><td rowspan="2">14.8 ± 0.4</td><td rowspan="2">0.11 ± 0.01</td><td rowspan="2">2050<sup>+35</sup><sub>-20</sub></td></tr><tr><td>2022</td><td>1556.3 ± 2.8</td><td>119.42 ± 0.11</td><td>8.99 ± 0.12</td><td>8.49 ± 0.13</td><td>15.6 ± 1.64</td><td>-107.4 ± 0.47</td></tr></tbody></table>

Notes. The mass ratio with respect to C is indicated as $q_{\text{Cb}}$.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Motion of Cb around C in RA/Dec for the 1000 best draws from our posterior distributions. We note that LLH stands for log likelihood.</figcaption></figure>

## 4. Discussion and concluding remarks

Before BEAST began, no $\le 30 M_{\odot}$ $M_{\text{Jup}}$ companion was known around stars more massive than $3 M_{\odot}$ – with only sporadic detections by RV in the mass range 2.5–3 $M_{\odot}$ (Quirrenbach et al. 2019; Wolthoff et al. 2022), questioning their very existence (see discussion in Janson et al. 2021a). The discovery of a circumstellar ~11 $M_{\text{Jup}}$ planet around the 6–10 $M_{\odot}$ binary b Centauri (Janson et al. 2021b) and one (possibly two) brown dwarfs close to the deuterium-burning limit around the ~9 $M_{\odot} \mu^2$ Scorpii (Squicciarini et al. 2022) first provided evidence for such a population, opening up a plethora of questions about its genesis.

The architecture of the HIP 81208 system turns out to be unique in many respects (Table 3). Not only is the B-type primary surrounded by a brown dwarf and an M-type stellar companion, but the additional discovery presented in this Letter of a ~15 $M_{\text{Jup}}$ companion to the C component makes it the first binary system with substellar companions to both components ever discovered by imaging.

Even if considered in isolation, a ≈15 $M_{\text{Jup}}$ companion at ~20 au from a late M-type star such as HIP 81208 C would be deemed remarkable. Figure 6 shows the mass ratio of confirmed giant planets and brown dwarfs ($M \in [1, 80] M_{\text{Jup}}$) around late M-type stars ($M_* \in [0.08, 0.3] M_{\odot}$): among five such DI companions, only two – WISE J0720-0846 (Burgasser et al. 2015) and LHS 2397a B (Dupuy et al. 2009) – have orbital separations <50 au; both of them, however, are characterized by a much larger mass ratio ($q \approx 0.7$) than $q_{\text{Cb}}$, which is indicative of a binary-like formation (Lodato et al. 2005). By including indirect techniques (with $a \in [3, 50]$ au), only two $a < 10$ au objects discovered via microlensing – OGLE-2016-BLG-0263L b (Han et al. 2017) and OGLE-2013-BLG-0911L b

<table><caption>**Table 3.** Emerging architecture of the hierarchical HIP 81208 system, with the median and (16th, 84th) percentiles used for *a* and *T*.</caption><thead><tr><th>Body</th><th>Primary</th><th><i>M</i> (<i>M</i><sub>☉</sub>)</th><th><i>a</i> (au)</th><th><i>T</i> (yr)</th></tr></thead><tbody><tr><td>A</td><td>A</td><td>2.58 ± 0.06</td><td>—</td><td>—</td></tr><tr><td>B</td><td>A</td><td>0.064<sup>+0.006</sup><sub>−0.007</sub></td><td>53.98<sup>+32.22</sup><sub>−15.00</sub></td><td>246.9<sup>+251.3</sup><sub>−95.4</sub></td></tr><tr><td>C</td><td>A</td><td>0.135<sup>+0.010</sup><sub>−0.013</sub></td><td>234.27<sup>+168.65</sup><sub>−68.96</sub></td><td>2232.4<sup>+4429.4</sup><sub>−1213.6</sub></td></tr><tr><td>Cb</td><td>C</td><td>0.0141 ± 0.0004</td><td>23.04<sup>+13.88 (a)</sup><sub>−6.55</sub></td><td>285.00<sup>+293.67</sup><sub>−112.07</sub></td></tr></tbody></table>

Notes. Values taken from Viswanath et al. (2023) for B and C. $^{(a)}$We note that the distribution of the semi-major axis of Cb is not Gaussian and the peak of the distribution is at about 17 au (see text). The column “primary” indicates the parent body of each object, to which *a* and *T* also refer.

(Miyazaki et al. 2020), both with a small $q \approx 0.03$ – were added$^4$.

While a full formation analysis of HIP 81208 Cb is beyond the scope of this work, it is worth mentioning possible formation pathways for the object and the whole quadruple system. Pivotal to a full understanding of the observed architecture is the formation of HIP 81208 C: the M star could be an outcome of either turbulent fragmentation of a star-forming core (Offner et al. 2010) – possibly followed by inward migration (Kuffmeier et al. 2019) – or GI within the disk of HIP 81208 A (Kratter et al. 2010), with the rough cutoff between these mechanisms (~500 au) likely depending on the environment and stellar mass (Offner et al. 2016).

If the former is true, the circumstellar disks of A and C would be truncated due to mutual gravitational actions (Panić et al. 2021). Provided no significant alteration of initial orbital parameters, a tentative estimate of the truncation radii $R_T$ for the two stars could be derived as in Artymowicz & Lubow (1994) by drawing $10^6$ ($a, e$) values for the A-C pair from the corresponding posterior distributions:

$$R_{T,C} = 0.88 \cdot R_R = 0.88 \cdot r_p \cdot \frac{0.49 \cdot q^{2/3}}{0.6 \cdot q^{2/3} + \ln(1 + q^{1/3})} = 26_{-9}^{+16} \text{ au}, \quad (2)$$

and

$$R_{T,A} = 0.88 \cdot (r_p - R_R) = 130_{-40}^{+70} \text{ au}, \quad (3)$$

where $r_p$ indicates the periastron of the orbit, and the Roche lobe $R_R$ is derived from the empirical formula by Eggleton (1983). The current position of Cb would be only marginally within the locations of its alleged parent disk, whence it would have sprouted either via CA (Mordasini et al. 2012; Alibert et al. 2013; Emsenhuber et al. 2021) or GI. According to CA models, the formation of a high-$q$ 15 $M_{\text{Jup}}$ object around a late M-type star is not expected (Kennedy & Kenyon 2008;

<footnote>$^4$ Data from the Exoplanet Encyclopaedia (http://exoplanet.eu/).</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Mass ratio vs. stellar mass for confirmed giant planets and brown dwarf companions ($1 M_{\text{Jup}} < M < 80 M_{\text{Jup}}$) to late M-type and massive brown dwarf primaries ($0.03 M_{\odot} < M < 0.3 M_{\odot}$). Masses of RV companions actually correspond to $M \sin i$. Detection methods are labeled with different colors; the semimajor axes or, if not available, projected separations have been grouped in four bins, with each one being associated with a circle size. The position of HIP 81208, provided with error bars, is explicitly labeled.</figcaption></figure>

Adams et al. 2021; Schlecker et al. 2022), which is also due to formation timescales (~$10^{6-7}$ Myr) exceeding typical disk lifetimes by one or two order of magnitude at separations ≥10 au (Dodson-Robinson et al. 2009). Conversely, GI could represent a formation channel for Cb (Meru & Bate 2010; Forgan & Rice 2013; Kratter & Lodato 2016), provided an unusually large disk-to-star ratio of ~30% ($\approx 40 M_{\text{Jup}}$) compared to the expected 1–10% (Manara et al. 2018; Haworth et al. 2020; Mercer & Stamatellos 2020, and references therein). Interestingly enough, the C-Cb separation would be within the typical range of M-type stellar binaries (Winters et al. 2019; Susemichel & Meyer 2022).

According to the alternative scenario, C itself would have formed via GI within the disk of A: simulations are able to produce objects with masses as high as $0.12 M_{\odot}$ already around $<1.2 M_{\odot}$ stars (Forgan et al. 2018). Despite the lack of GI models suited to B-type hosts, evidence for a more-than-linear dependence of disk mass on stellar mass (Pascucci et al. 2016), coupled to observed circumstellar disks to late B stars spanning hundreds of au (Alonso-Albi et al. 2009; Mariñas et al. 2011; Garufi et al. 2017), tentatively hints toward such a possibility. As a companion to a disk-borne object, HIP 81208 Cb would intriguingly retain – whether it formed in situ or via dynamical capture (Podsiadlowski et al. 2010; Ochiai et al. 2014) – the hierarchical level of a satellite (Lazzoni et al. 2022).

A detailed characterization of the orbital parameters of B, C, and Cb, and in particular their mutual inclinations, will discriminate between the two scenarios, helping shed light on this unique multiple system.

**Acknowledgements.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (COBREX; grant agreement n° 885593). SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zürich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European

Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 Exploring the Diversity of Extrasolar Planets (MA 9185/1) and from the Swiss National Science Foundation under grant 200021_204847 PlanetsInTime. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. R.G. and S.D. acknowledge the support of PRIN-INAF 2019 Planetary Systems At Early Ages (PLATEA). This research has made use of the SIMBAD database and VizieR catalogue access tool, operated at CDS, Strasbourg, France. This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project of the Univ. Grenoble-Alpes. This work was supported by the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES. This research has made use of data obtained from or tools provided by the portal http://exoplanet.eu/ of The Extrasolar Planets Encyclopaedia.

## References

* Adams, F. C., Meyer, M. R., & Adams, A. D. 2021, ApJ, 909, 1

* Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109

* Allard, F. 2016, in SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, J. Richard, L. Cambrésy, et al., 223

* Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357

* Alonso-Albi, T., Fuente, A., Bachiller, R., et al. 2009, A&A, 497, 117

* Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651

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

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

* Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391

* Boss, A. P. 1997, Science, 276, 1836

* Bouwman, J., Lawson, W. A., Dominik, C., et al. 2006, ApJ, 653, L57

* Brandeker, A., & Cataldi, G. 2019, A&A, 621, A86

* Burgasser, A. J., Melis, C., Todd, J., et al. 2015, AJ, 150, 180

* Cantalloube, F., Dohlen, K., Milli, J., Brandner, W., & Vigan, A. 2019, The Messenger, 176, 25

* Chauvin, G., Videla, M., Beust, H., et al. 2023, A&A, 675, A114

* Chomez, A., Lagrange, A.-M., Delorme, P., et al. 2023, A&A, 675, A205

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

* Cuntz, M., Luke, G. E., Millard, M. J., Boyle, L., & Patel, S. D. 2022, ApJS, 263, 33

* Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246

* Damiani, F., Prisinzano, L., Pillitteri, I., Micela, G., & Sciortino, S. 2019, A&A, 623, A112

* Deeg, H. J., & Belmonte, J. A. 2018, Handbook of Exoplanets (Springer)

* Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al.

* Desidera, S., Bonomo, A. S., Claudi, R. U., et al. 2014, A&A, 567, L6

* de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354

* Dodson-Robinson, S. E., Veras, D., Ford, E. B., & Beichman, C. A. 2009, ApJ, 707, 79

* Dohlen, K., Langlois, M., Saisse, M., et al. 2008, SPIE Conf. Ser., 7014, 70143L

* Dupuy, T. J., Liu, M. C., & Ireland, M. J. 2009, ApJ, 699, 168

* Eggleton, P. P. 1983, ApJ, 268, 368

* Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A70

* Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138

* Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020a, A&A, 634, A2

* Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020b, A&A, 637, A9

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

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Forgan, D., & Rice, K. 2013, MNRAS, 432, 3168

* Forgan, D. H., Hall, C., Meru, F., & Rice, W. K. M. 2018, MNRAS, 474, 5036

* Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1
Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92

Garufi, A., Meeus, G., Benisty, M., et al. 2017, A&A, 603, A21

Han, C., Udalski, A., Gould, A., et al. 2017, AJ, 154, 133

Hartkopf, W. I., & Mason, B. D. 2004, Rev. Mex. Astron. Astrophys. Conf. Ser., 2, 83

Haworth, T. J., Cadman, J., Meru, F., et al. 2020, MNRAS, 494, 4130

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

Janson, M., Gratton, R., Rodet, L., et al. 2021a, Nature, 600, 231

Janson, M., Squicciarini, V., Delorme, P., et al. 2021b, A&A, 646, A164

Kennedy, G. M., & Kenyon, S. J. 2008, ApJ, 673, 502

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585

Kuffmeier, M., Calcutt, H., & Kristensen, L. E. 2019, A&A, 628, A112

Lazzoni, C., Desidera, S., Gratton, R., et al. 2022, MNRAS, 516, 391

Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44

Lodato, G., Delgado-Donate, E., & Clarke, C. J. 2005, MNRAS, 364, L91

Lueber, A., Kitzmann, D., Fisher, C. E., et al. 2023, ApJ, submitted [arXiv:2305.07719]

Maire, A.-L., Langlois, M., Delorme, P., et al. 2021, J. Astron. Telesc. Instrum. Syst., 7, 035004

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

Mariñas, N., Telesco, C. M., Fisher, R. S., & Packham, C. 2011, ApJ, 737, 57

Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Mercer, A., & Stamatellos, D. 2020, A&A, 633, A116

Meru, F., & Bate, M. R. 2010, MNRAS, 406, 2279

Miret-Roig, N., Bouy, H., Raymond, S. N., et al. 2022, Nat. Astron., 6, 89

Miyazaki, S., Sumi, T., Bennett, D. P., et al. 2020, AJ, 159, 76

Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111

Musielak, Z. E., Cuntz, M., Marshall, E. A., & Stuit, T. D. 2005, A&A, 434, 355

Ochiai, H., Nagasawa, M., & Ida, S. 2014, ApJ, 790, 92

Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 725, 1485

Offner, S. S. R., Dunham, M. M., Lee, K. I., Arce, H. G., & Fielding, D. B. 2016, ApJ, 827, L11

Panić, O., Haworth, T. J., Petr-Gotzens, M. G., et al. 2021, MNRAS, 501, 4317

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59

Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38

Podsiadlowski, P., Rappaport, S., Fregeau, J. M., & Mardling, R. A. 2010, ArXiv e-prints [arXiv:1007.1418]

Quirrenbach, A., Trifonov, T., Lee, M. H., & Reffert, S. 2019, A&A, 624, A18

Ratzenböck, S., Großschedl, J. E., Alves, J., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202346901

Roberts, Lewis C., J., Tokovinin, A., Mason, B. D., et al. 2015, AJ, 149, 118

Schlecker, M., Burn, R., Sabotta, S., et al. 2022, A&A, 664, A180

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

Squicciarini, V., & Bonavita, M. 2022, A&A, 666, A15

Squicciarini, V., Gratton, R., Bonavita, M., & Mesa, D. 2021, MNRAS, 507, 1381

Squicciarini, V., Gratton, R., Janson, M., et al. 2022, A&A, 664, A9

Susemiehl, N., & Meyer, M. R. 2022, A&A, 657, A48

Udry, S., Dumusque, X., Lovis, C., et al. 2019, A&A, 622, A37

Viswanath, G., Janson, M., Gratton, R., et al. 2023, A&A, 675, A54

Wang, J., Fischer, D. A., Xie, J.-W., & Ciardi, D. R. 2014, ApJ, 791, 111

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

Wolthoff, V., Reffert, S., Quirrenbach, A., et al. 2022, A&A, 661, A63

Yuan, H. B., Liu, X. W., & Xiang, M. S. 2013, MNRAS, 430, 2188

Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2003, A&A, 404, 775

Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2004, A&A, 426, 695

<footnote n="1">¹ LESIA, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France
e-mail: antoine.chomez@obspm.fr</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS-INSU, Institut de Planetologie et d'Astrophysique de Grenoble (IPAG), UMR 5274, Grenoble 38041, France</footnote>

<footnote n="3">³ INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="4">⁴ Institutionen för astronomi, Stockholms universitet; AlbaNova universitetscentrum, 106 91 Stockholm, Sweden</footnote>

<footnote n="5">5 Univ. Lyon, Univ. Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon (CRAL) UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="6">6 Université Côte d'Azur, OCA, CNRS, Lagrange, 96 Bd de l'Observatoire 06300 Nice, France</footnote>

<footnote n="7">7 Pixyl S.A., 5 Av. du Grand Sablon, 38700 La Tronche, France</footnote>

<footnote n="8">8 ETH Zurich, Institute for Particle Physics and Astrophysics, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote n="9">9 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="10">10 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="11">11 Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany</footnote>

<footnote n="12">12 Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany</footnote>

<footnote n="13">13 Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, 3012 Bern, Switzerland</footnote>

<footnote n="14">14 Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109, USA</footnote>

<footnote n="15">15 Landessternwarte, Zentrum für Astronomie der Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany</footnote>

<footnote n="16">16 Université Grenoble Alpes, CNRS, Observatoire des Sciences de l'Université de Grenoble (OSUG), 122 Rue de la Piscine, 38400 Saint-Martin-d'Hères, France</footnote>
**Appendix A: Details on the custom PSF subtraction**
**routine**

The custom PSF subtraction routine, sketched in Sect 2.2, pro-
ceeded in two steps based on small stamps of 47×47 pixels
roughly centered on C. In the first step, we used a 47x47 pixel
stamp of the mean of the two off-axis PSFs as a model, acquired
just before and just after the observations. We recursively removed
C and then Cb from each of the 48 individual frames of the coron-
agraphic sequence by minimizing the residuals between the data
stamp and our empirical PSF model stamp. This minimization was
performed by injecting negative versions of this empirical PSF
model in a local grid centered on C oversampled by a factor 100
in each spatial dimension. We used three free parameters each for
C and Cb, namely the oversampled pixels' positions in x and y
and the source to model contrast, and we selected the models that
minimize the absolute value of the residuals on small optimization
zones (1000 by 1000 oversampled pixels). After this first step, we
noticed that the residuals after subtraction of this first PSF model
were characterized by the following: 1) a relatively high intensity,
on the order of 1-2% of the local flux of C, and 2) a systematic
shape independent of time but with an alignment following the
parallactic angle rotation. We interpreted these features as hints
that the local PSF at C's position slightly differed from the cali-
bration off-axis PSF. We therefore added C again (using its fitted
parameters) on each of the 48 stamps with both C and Cb removed,
effectively building a local PSF with the same pupil orientation on
each frame. Afterwards, we applied a similar approach to stan-
dard ADI, which median-combined the 48 resulting subframes
obtained at different parallactic angles without derotating them,
producing a local pupil-stabilized PSF. As in ADI, the weak resid-
uals from the subtraction of Cb, already close to the background
noise and rotating around C with the parallactic angle, were fur-
ther removed from this local PSF model by means of the nondero-
tated median.

As a second step, we repeated the same minimization
approach, starting from the 48 small stamps containing both C
and Cb, but using this local pupil-stabilized PSF as a model
instead of the off-axis PSF. Minimization directly provides the

best-fit parameters for position and contrast for C and Cb on
each frame, albeit in a pupil-tracking rotating frame of reference.
These measurements were then derotated to sky coordinates and
averaged to obtain the results shown in Table C.1. The intensity
of the residuals was reduced by a factor ~5 after this second step
compared to the first step, with the absolute value of the highest
residuals in any individual frames at both epochs being barely
above the local background noise.

We determined the uncertainties associated with our esti-
mates of positions and flux by deriving the standard devia-
tion of the resulting 48 independent measurements, which nat-
urally and robustly include all sources of systematic errors that
cause frame to frame variations, such as tip-tilt jitter or atmo-
spheric transmission variability during the observing sequence.
The main remaining systematic – namely the uncertainty on the
position of the central star, which is constant over the sequence
– was quadratically added to the measured uncertainties and
dominates the error budget. In the peculiar case of the rela-
tive position of Cb around C (reported in Table 2), this system-
atic is naturally canceled out and the dynamical fits performed
in subsection 3.3 consequently employ the much smaller
error bars obtained when this systematic contribution is taken
out.

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Residuals after subtracting both the PSF of C and Cb (with 2019 on the left and 2022 on the right) for the K1 band. The position of C (subtracted) is highlighted by the white cross. Residuals are below 1 ADU.</figcaption></figure>
**Appendix B: Raw coronagraphic frames displaying**
**the ADI rotation around C of Cb**

As mentioned in Sect. 2.3, the presence of the source Cb is visu-
ally evident even in raw coronagraphic frames despite its close

vicinity to HIP 82108 C. A subset of the 48 individual frames is
provided in Fig. B.1.

<figure><img src="image_9.png" /><figcaption>**Fig. B.1.** Collection of raw coronagraphic frames of the 2019 dataset. The current position of the Cb companion at the given frame is highlighted by the white arrow, and past positions (projected onto the current frame) are highlighted by the white crosses. We can clearly see the expected rotation due to the ADI sequence. A PSF artifact would have resided at a fixed angle relative to C in all frames. (a) First raw coronagraphic frame of the 2019 dataset (b) 25th raw coronagraphic frame of the 2019 dataset (c) Last raw coronagraphic frame of the 2019 dataset</figcaption></figure>
**Appendix C: Astrophotometric results**

We provide in Table C.1 a summary of the astrometric and photometric results obtained for HIP 81208 Cb (Sec. 2.3) by means of

the PACO and PSF subtraction reductions described in Sec. 2.2. Values for the newly found background source CC14 – only detected by PACO – are provided in Table C.2.

<table><caption>**Table C.1.** Astrometry and photometry extracted for the companions B, C, and Cb with PACO and the PSF subtraction.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="3">2019-08-06</th><th colspan="3">2022-04-05</th></tr><tr><th>HIP 81208 B</th><th>HIP 81208 C</th><th>HIP 81208 Cb</th><th>HIP 81208 B</th><th>HIP 81208 C</th><th>HIP 81208 Cb</th></tr></thead><tbody><tr><td rowspan="4">PACO ASDI</td><td>Sep (mas)</td><td>315.5±3.1</td><td>1495.1±3.2</td><td>1573.4±3.2</td><td>324.4±2.7</td><td>1493.0±2.8</td><td>1556.3±2.8</td></tr><tr><td>PA (°)</td><td>357.6±0.57</td><td>116.3±0.13</td><td>119.9±0.12</td><td>1.39±0.48</td><td>115.95±0.09</td><td>119.42±0.11</td></tr><tr><td>ΔK1 (mag)</td><td>6.99±0.07</td><td>5.85±0.08</td><td>8.93±0.08</td><td>6.91±0.11</td><td>5.93±0.12</td><td>8.99±0.12</td></tr><tr><td>ΔK2 (mag)</td><td>6.70±0.08</td><td>5.48±0.08</td><td>8.59±0.08</td><td>6.60±0.16</td><td>5.49±0.12</td><td>8.49±0.13</td></tr><tr><td rowspan="4">PSF subtraction</td><td>Sep (mas)</td><td>-</td><td>1494.4±3.2</td><td>1573.1±3.3</td><td>-</td><td>1491.4±2.8</td><td>1554.8±3.2</td></tr><tr><td>PA (°)</td><td>-</td><td>116.36±0.14</td><td>119.76±0.15</td><td>-</td><td>115.89±0.15</td><td>119.2±0.16</td></tr><tr><td>ΔK1 (mag)</td><td>-</td><td>5.80±0.03</td><td>8.96±0.04</td><td>-</td><td>5.81±0.03</td><td>8.97±0.08</td></tr><tr><td>ΔK2 (mag)</td><td>-</td><td>5.54±0.03</td><td>8.55±0.04</td><td>-</td><td>5.48±0.03</td><td>8.45±0.08</td></tr></tbody></table>

**Notes:** The astrometric values are averaged over the two wavelengths.

<table><caption>**Table C.2.** Astrometry and photometry for the additional background source CC14 detected in this analysis.</caption><thead><tr><th></th><th>Sep (mas)</th><th>PA (°)</th><th>ΔK1</th><th>ΔK2</th></tr></thead><tbody><tr><td>2019-08-06</td><td>3923.3±6.7</td><td>141.1±0.1</td><td>14.57±0.26</td><td>14.03±0.27</td></tr><tr><td>2022-04-05</td><td>3875.5±8.14</td><td>140.4±0.13</td><td>15.87±0.28</td><td>15.04±0.26</td></tr></tbody></table>
## Appendix D: The bound nature of A-B and C-Cb

Whilst the proper motion analysis of Cb firmly allowed us to conclude on its nonbackground nature and its common motion to A, B, and C, it does not exclude, in principle, an alternative hypothesis:

The A-B system is totally independent of the C-Cb system, both being Sco-Cen binaries projected by chance at a short separation from one another.

In order to quantify the probability of this alternative scenario, we adapted the argument already made for HIP 81208 B (Viswanath et al. 2023) and $\mu^2$ Scorpii b (Squicciarini et al. 2022). After defining indicative coordinate limits for UCL as $(l, b) = [313^\circ, 343^\circ] \times [2^\circ, 28^\circ]$, we recovered $N_0 = 3835$ bona fide members to this subgroup from the Gaia DR2-based list of Sco-Cen members assembled by Damiani et al. (2019). At the distance and age of Sco-Cen, the census of the stellar population of the association is reasonably complete$^5$. However, a source can be overlooked by Gaia if it happens to be located too close to a brighter star, that is, if the $\Delta G$ between the former and the latter is larger than the maximum contrast achievable by the satellite at the corresponding angular separation $s$. We therefore define, as a shaded area the circular region centered on a star within which the average detection efficiency $\delta(s, \Delta G)$ of Gaia equals 50% for a given apparent G magnitude, and effective separation $s_{\text{eff}}$.

Our goal here is to quantify the number of these "phantom" UCL stars, so as to enable an estimation of the probability of spotting at least one of them within the entire BEAST survey. Intuitively, the computation hinges upon the following: 1) the total shaded area $A_s$, obtained as the sum of individual shaded areas for all Gaia sources within the boundaries of UCL; and 2) the number – corrected for completeness – of UCL members, $N_{\text{UCL}}$. In regards to the former, we queried Gaia DR3 (Gaia Collaboration 2023), finding approximately $8 \cdot 10^7$ stars within the coordinate limits of UCL. From Brandeker & Cataldi (2019), we then recovered the detection efficiency of Gaia DR2 as a function of $\Delta G$ and $s$, $\delta(s, \Delta G)$. In this way, we were able to compute, for every Gaia source $i$, the effective separation $s_{\text{eff}}$ as a function of the apparent G magnitude of a hypothetical phantom star, $G$:

$$s_{\text{eff},i} = s | \bar{\delta}(s, \Delta G_i) = 0.5, \qquad (\text{D.1})$$

where

$$\bar{\delta}(s, \Delta G_i) = \frac{1}{s^2} \int_0^s \delta(\bar{s}, \Delta G_i) \bar{s}^2 d\bar{s} \qquad (\text{D.2})$$

and $\Delta G_i(G) = G - G_i$. Summation over Gaia stars yields the total shaded area as a function of $G$:

$$A_s(G) = \sum_i \pi s_{\text{eff},i}^2. \qquad (\text{D.3})$$

The probability density function (PDF) of phantom stars can be now expressed as follows:

$$n_P(G) = \frac{A_s(G)}{A_{\text{UCL}} - A_s(G)} \cdot N_{\text{UCL}} \cdot 0.5 \cdot \zeta_{\text{Gaia}}(G), \qquad (\text{D.4})$$

where $\zeta_{Gaia}(G)$ is the apparent G magnitude PDF of UCL stars, and the factor 0.5 accounts for the expectation that 50% of these object were already detected by Gaia.

To cope with the incompleteness of the initial mass function (IMF) of the UCL sample at its faint end (i.e., for unseen substellar objects), we recovered the sample of Upper Scorpius$^6$ (US) members by Miret-Roig et al. (2022), pushing completeness down to ~10$M_{\text{Jup}}$. 2MASS J magnitudes were converted into Gaia G magnitudes based again on BHAC15 isochrones at 15 Myr, yielding the $\zeta_{US}(G)$ PDF of the sample; a new normalized $\zeta(G)$ could then be built by combining the Gaia-based UCL list and the US sample, setting a sharp transition between the former and the latter at the dimmest magnitude where they intersect ($\hat{G} = 16.8$ mag). Above this value, the two distributions start to significantly differ due to Gaia incompleteness. The number of unseen sources $n_U$ recovered in this way amounts to $\approx 700$ (~18%), and their PDF is given by the following:

$$n_U(G) = \begin{cases} 0 & G < \hat{G} \\ N_{\text{UCL}} \cdot [\zeta_{\text{US}}(G) - \zeta_{\text{Gaia}}(G)] & G \ge \hat{G} \end{cases} \qquad (\text{D.5})$$

In order to compute the probability that a phantom or an unseen star is hidden by a BEAST star, we consider as a typical BEAST star an object as bright in the apparent G band as the mean of the sample ($G = 5.29$); for a contrast $\Delta G \approx 8$ mag, the effective shading separation $s_{\text{eff,B}}$ of this star starts being larger than the half-edge (5.5") of IRDIS FOV ($A_{\text{IRD}} = 11'' \times 11''$); we therefore impose $s_{\text{eff,B}}(\Delta G) = \inf(s_{\text{eff,B}}(\Delta G), \sqrt{A_{\text{IRD}}/\pi})$. The differential probability associated with the event as a function of $G$ is given by

$$f(G) = f_P(G) + f_U(G) = \frac{n_P(G) \cdot \pi s_{\text{eff,B}}^2(G)}{A_s(G)} + \frac{n_U(G) \cdot A_{\text{IRD}}}{A_{\text{UCL}}}, \qquad (\text{D.6})$$

where the second term takes into account that unseen sources should be spread over the entire UCL. Integration of $f(G)$ yields $p = \int_{G=5.29}^{G=25} f(G)dG = 2.8 \cdot 10^{-5}$. The false alarm probability associated with the event of finding at least one such object across the whole survey, having completed the observations of 47 stars as of yet, is equal to the following:

$$\text{FAP} = 1 - \binom{47}{0} (1-p)^{47} = 1.3 \times 10^{-3}. \qquad (\text{D.7})$$

In order to evaluate the impact of the assumption of a constant age for US, which is instead known to have experienced a long-lasting star formation history ranging between 15 and 5 Myr ago (Squicciarini et al. 2021; Ratzenböck et al. 2023), the conversion of J magnitudes into G magnitudes was repeated by supposing a constant age 5 Myr. The resulting FAP $\approx 1.4 \cdot 10^{-3}$ provides robust evidence for the independence of the result on model assumptions, firmly allowing us to exclude the alternative scenario in favor of the one positing a single quadruple system.

<footnote>$^5$ According to BHAC15 isochrones (Baraffe et al. 2015) at solar metallicity, a 0.08 $M_\odot$ star aged 15 Myr has absolute $G = 11.45$, corresponding to an apparent $G \sim 17$ mag at the mean separation of UCL (~140 pc); the survey is virtually complete for $G \in [12, 17]$ mag (Gaia Collaboration 2018).</footnote>

<footnote>$^6$ Upper Scorpius is, together with UCL and Lower Centaurus-Crux, one of the three subgroups in which Sco-Cen is classically divided (de Zeeuw et al. 1999). We verified through a Kolmogorov-Smirnov test with $\alpha = 0.05$ that the absolute G-magnitude distribution for US stars – selected from the same sample adopting coordinate boundaries as in Squicciarini et al. (2021) – is consistent with its UCL analog for $G > 4$.</footnote>
# Appendix E: Characterization of Cb

The derivation of a photometric mass estimate for HIP 81208 Cb is mediated by MADYS (Squicciarini & Bonavita 2022), as in previous BEAST publications (Squicciarini et al. 2022; Viswanath et al. 2023). After averaging $K1$ and $K2$ contrasts derived by PACO over the two epochs, the conversion of those contrasts into calibrated apparent magnitudes was operated by means of the 2MASS $K_s$ magnitude of the primary. Being HIP 81208 A classified as a B9V star, the impact of the approximation $K_{s,A} \approx K1_A \approx K2_A$ is well within the photometric error budget (Pecaut & Mamajek 2013).

Interstellar reddening toward HIP 81208 is known to be rather small (Viswanath et al. 2023) and translating, in the $K_s$ band, to a negligible $A_{K_s} = 0.003 \pm 0.006$ adopting $A_{K_s}/E(B-V) = 0.306$ (Yuan et al. 2013). Likewise, the adopted parallax and age estimates reflect those used in the paper by Viswanath et al. (2023).

Having obtained absolute magnitudes for the object,

$$K1 = 9.90 \pm 0.04 \text{ mag},$$

$$K2 = 9.48 \pm 0.04 \text{ mag},$$

we built a set $M$ of substellar evolutionary models that are adequate for the age and mass range of interest, while providing synthetic SPHERE magnitudes at once: such a set, to which observed magnitudes were compared, includes the Ames-Dusty models (Allard et al. 2001) and the BT-Settl models (Allard 2016). Details on the fitting algorithm, encompassing all sources of uncertainty within a Monte Carlo framework, can be found in Squicciarini & Bonavita (2022). The output of each model $i \in M$ corresponds to a triplet $(M_{min,i}, M_{opt,i}, M_{max,i})$ equal to the $(16^{th}, 50^{th}, 84^{th})$ percentiles of the posterior mass distribution; the two outputs were averaged in the following manner:

$$M_{min} = \inf_{i \in M} (M_{min,i}) \qquad (E.1)$$

$$M_{opt} = \text{mean}(\{M_{opt,i}\}) \qquad (E.2)$$

$$M_{max} = \sup_{i \in M} (M_{+,i}) \qquad (E.3)$$

with the goal of embedding theoretical uncertainties onto the final estimate.

The posterior mass distribution returned by each model can be easily translated into the posterior distribution of any astrophysical parameter of interest provided by the isochrone grids. We were therefore able to derive in a self-consistent way (using similar equations to Eq. E.1-E.3) the best-fit estimates for effective temperature, surface gravity, and bolometric luminosity. Likewise, we also computed synthetic 2MASS $H$ and $K$ magnitudes as a helpful first-guess estimate for follow-up studies:

$$H = 10.28 \pm 0.07 \text{ mag},$$

$$K_s = 9.74 \pm 0.04 \text{ mag}.$$

We highlight that independent mass determinations were derived a posteriori, and used as a control sample, starting from the best-fit $\log L/L_{\odot}$ through the recent ATMO2020 (Phillips et al. 2020) and Sonora Bobcat (Marley et al. 2021) grids$^7$; these best-fit masses of $14.51_{0.15}^{+0.16} M_{\text{Jup}}$ and $14.44_{0.14}^{+0.15} M_{\text{Jup}}$, respectively, are consistent with our best-fit mass estimates. Nonetheless, as already mentioned in Sect. 3.2, we are not able to exclude the possibility of unaccounted systematic effects, which are common to all the adopted models.

# Appendix F: Orbital fit of Cb: Corner plot

Given the small separation between C and Cb, most sources of systematic error are either canceled out (centering error) or significantly decreased (platescale and true north error), enabling an accurate determination of their relative separation at both epochs (Appendix A).

Starting from the relative C-Cb astrometry as measured by PSF subtraction and their best-fit masses (cp. Table 3), we ran an MCMC code based on the emcee code (Foreman-Mackey et al. 2013) in order to derive the orbital parameters of Cb's orbit around C.

The input parameters for the MCMC are a logarithmically uniform prior for the semimajor axis ($a \in [0, 80]$ au) and an eccentricity $e \in [0, 0.4]$.

The posterior distribution for the orbital parameters of Cb derived in this work is provided in Fig. F.1.

<footnote>$^7$ These grids were not included in $M$ because they are currently not equipped with SPHERE filters.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. F.1.** Corner plot of the orbital solution for Cb. Reading from left to right, the acronyms stand for (from top to bottom) semimajor axis, eccentricity, argument of periapsis, periapsis time, orbital period, Cb mass, inclination, and longitude of the ascending node.</figcaption></figure>