© The Authors 2023

1 Introduction

Most known exoplanets are detected through transit and radial velocity techniques because of the exquisite sensitivity of these latter to planets in small orbits around their host stars (e.g. Mayor & Queloz 1995; Winn et al. 2011). However, both techniques are much less sensitive to planets in wide orbits (several au or larger), and therefore this class of planets has been much less explored. In recent years, a growing population of wide planets has been discovered with direct imaging (e.g. Macintosh et al. 2015; Chauvin et al. 2017; Keppler et al. 2018; Bohn et al. 2020) and, to an increasing extent, in synergy with astrometry (e.g. Brandt et al. 2021; Lacour et al. 2021; Hinkley et al. 2023; Franson et al. 2023; Currie et al. 2023). These techniques are also sensitive to brown dwarf companions (e.g. Janson et al. 2012; Crepp et al. 2016; Currie et al. 2020); while such companions are rare, their high masses and brightnesses make them suitable for detailed characterisation, which can be used to form important links in our understanding of atmospheric conditions down to progressively smaller substellar masses and temperatures.

Direct imaging surveys (e.g. Uyama et al. 2017; Nielsen et al. 2019; Launhardt et al. 2020; Vigan et al. 2021) have studied wide planet and brown dwarf companion demographics in a range of stellar system types, including binaries (e.g. Bonavita et al. 2016; Asensio-Torres et al. 2018), debris disk systems (e.g. Janson et al. 2013; Meshkat et al. 2017), and stars of different masses (e.g. Bowler et al. 2012; Delorme et al. 2012; Wagner et al. 2022). Studying substellar companion demographics as a function of stellar mass can yield insights into the planet and brown dwarf formation processes, and in this context, the most massive stars in the solar neighbourhood, namely B stars, represent an important part of the puzzle. To probe this regime deeply for the first time, we initiated the B-star Exoplanet Abundance Study (BEAST; see Janson et al. 2021a), which is a high-contrast imaging survey using the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2019) Extreme Adaptive Optics instrument at the Very Large Telescope (VLT) in Chile to study 85 B stars in the nearby young Sco-Cen region (de Zeeuw et al. 1999). Currently ongoing, BEAST has already discovered several new planets and low-mass substellar companions (Janson et al. 2019, 2021b; Squicciarini et al. 2022).

Here, we report the discovery of two low-mass companions, both in a dynamically unusual orbital configuration around the B9 star HIP 81208. The paper is structured as follows: in Sect. 2, we give the details of the observation of the target and subsequent data reduction. In Sect. 3, we describe the target star and derive some of its basic parameters. Section 4 describes the main results from this work, including characterisation of the candidate properties, a detailed analysis of the spectrum of the inner companion candidate, and analysis of the orbits of the two companion candidates, as well as a discussion on the dynamical stability of the system. The main conclusions from this study are outlined in Sect. 5.

2 Observations and data reduction

The target system, HIP 81208, was observed twice with the SPHERE instrument at the Very Large Telecsope (VLT), both in the context of BEAST. The first observation was acquired on 6 Aug 2019 (MJD 58701.02), and a second-epoch follow-up observation was acquired on 5 April 2022 (MJD 59674.35). Both observations used the standard settings for BEAST (Janson et al. 2021a) in the IRDIFS-EXT mode (Zurlo et al. 2014), which allows simultaneous low-resolution YJH-band spectroscopy using the Integral Field Spectrograph (IFS; Beuzit et al. 2019) and dual-band imaging in the K band with the K12 filter-pair (Vigan et al. 2010) using the Infra-Red Dual-band Imager and Spectrograph (IRDIS; Beuzit et al. 2019).

The main observational sequence was executed in pupil-stabilised mode – which facilitates angular differential imaging (ADI, see Marois et al. 2006) reductions – and utilised the N-ALC-YJH-S coronagraph to enhance the instrumental contrast close to the parent star. This long coronagraphic sequence was complemented by two short pairs of images: A short unsaturated (non-coronagraphic) exposure of the star was acquired before and after the sequence for spectrophotometric calibration purposes, and a short coronagraphic exposure was taken with the so-called waffle mode (Cantalloube et al. 2019) turned on, which generates satellite images of the primary star for astrometric calibration purposes. An empty sky frame was also included in the observational procedure for the purpose of sky subtraction. The long coronagraphic sequence consisted of 16 frames, with three subintegrations of 64 s each. The total integration time on-target was therefore 3072 s in each epoch. The total field rotations during the observations were 51.9 deg in the 2019 epoch and 57.5 deg in the 2022 epoch.

Data reduction was performed according to the regular BEAST procedure (Janson et al. 2021a), with the SPHERE Data Center¹ (Delorme et al. 2017) software based on the SPHERE pipeline (Pavlov et al. 2008) for all basic calibration steps including dark and sky subtraction, flat fielding, and spectral extraction. For pixel scale, true north orientation and pupil offset, we used the long-term values given in Maire et al. (2021), which are based on a 5-yr analysis of the astrometric stability of SPHERE that relies on observations of star clusters taken at regular intervals. Accordingly, for both epochs, the pixel scale adopted for IRDIS was 12.258±0.004 mas pixel⁻¹ for K1 band and 12.253±0.003 mas pixel⁻¹ in the K2 band, while the pixel scale for IFS was 7.46 ± 0.02 mas pixel⁻¹. The true north orientation used was −1.77±0.04 deg and the pupil offset was 136± 0.03 deg. For ADI reduction purposes, we used SpeCal (Galicher et al. 2018) for IRDIS in Template Locally Optimised Combination of Images (TLOCI; Marois et al. 2014) mode, and a pipeline for IFS based on Karhunen-Loève Image Projection (KLIP, see e.g. Soummer et al. 2012).

3 Stellar properties

HIP 81208 (alias HD 149274, TYC 7357-207-1, TIC 280474618, Gaia DR3 6020514769906985728, 2MASS J16351384-3543287) is a young bright (V = 6.632 ± 0.006 mag) B9V star (Houk 1982) residing in the Upper Centaurus Lupus (UCL) subregion of the Sco-Cen stellar association (Hoogerwerf 2000), at a distance of ${148.7}_{-1.3}^{+1.5}\text{pc}$ from Earth. Some basic measurements of the star obtained from the literature are listed in Table 1.

The current literature gives no indication of a binary companion to HIP 81208. A potentially interesting Gaia star in the vicinity of HIP 81208 is discussed in Sect. 4.6. The proper motion anomaly (PMa) between the long-term HIPPARCOS-Gaia DR3 and short-term Gaia DR3 proper motion vectors is very low (PMaRA = 0.063 ± 0.056 mas yr⁻¹, PMaDE = 0.090 ± 0.041 mas yr⁻¹; see Kervella et al. 2022), ruling out the existence of any (>0.05 M_⊙) binary companion to the star within 3–30 au. There is also no hint of a faint or unresolved astro-metric binary companion within 2–20 au of the star in Gaia DR2, as inferred from the low value of the renormalised unit weight error (RUWE), of namely ρ = 0.898 (Belokurov et al. 2020). Additionally, Stock (2021) analysed the 26 XSHOOTER (Vernet et al. 2011) spectra available in the ESO archive and derived radial velocities with a precision of around 2 km s⁻¹. The radial velocities cover a period of slightly less than 8 yr from 2010 to 2017 and have an rms of about 5 km s⁻¹. The analysis revealed no indications of periodic variations, except for a slightly significant period of 1 yr, which can also be seen in the window function. From this, we conclude that there are no indications of spectroscopic companions to HIP 81208 in the XSHOOTER radial velocities. In the following subsections, we constrain the main parameters of the target, a summary of which is also listed in Table 2.

3.1 Extinction and reddening

Gaia DR3 lists a monochromatic extinction estimate of ${\text{A}}_0={0.0229}_{-0.0034}^{+0.0035}$ mag for the target at λ₀ = 541.4 nm (see Creevey et al. 2023, for the definition of A₀). A_V, the extinction in the Johnson V band centred at ~548 nm, is dependent on the spectrum of the emitting source in the V band and hence is intrinsically different from A₀, which is a property of the interstellar medium alone and is affected only by the amount of absorption therein. However, this dependence of A_V on the source spectrum is often negligible, especially for small extinctions (A₀ < 2 mag), at which A₀ and A_V follow an almost identity relation (see Fig. 11.4(a) in the Gaia DR3 online documentation²). Therefore, in this case, A₀ can be taken as a rough approximation of A_V for our target. The STILISM reddening maps³ from Lallement et al. (2018) at the position and distance of HIP 81208 gives a reddening estimate of E(B–V) = 0.011 ± 0.021 mag, which can be translated to an A_V value of 0.034 ± 0.064 mag using the relation A_V = 3.07×E(B–V) (McCall 2004, appropriate for lightly reddened A0 V stars). Examination of the reddening values for stars within ±4° of the position of HIP 81208 in the Reis et al. (2011) catalogue shows that stars in its vicinity are negligibly reddened until a distance of ~160 pc, and the stars in the distance range 170–180 pc have A_V = 0.1–0.2 mag. In particular, four stars (HD 147493, HD147387, HD147149, HD147597) in the immediate vicinity of HIP 81208, with distances of 146–150 pc, have negligible reddening values of E(B–V) = 0.000 ± 0.008, −0.008 ± 0.011, 0.003 ± 0.012, and 0.009 ± 0.011 mag, respectively, with an inferred A_V = 0.00 ± 0.03, −0.03 ± 0.05, 0.01 ± 0.05, and 0.04 ± 0.05 mag, respectively. Given the similar distances and positions of these stars in the sky to our target, one can expect them to have comparable reddening and extinction to HIP 81208.

An independent estimate of the reddening and extinction of this B star can be obtained from UBV photometry using the Q-method to calculate the reddening-free Q index (e.g. Johnson & Morgan 1953, see Appendix C.3 of Pecaut & Mamajek 2013). Adopting the U–B colour from Slawson et al. (1992) and the B–V colour from ESA (1997), and the updated calibration of the Q-index to the dereddened (B–V)₀ and (U–B)₀ colours from Pecaut & Mamajek (2013), we estimate the intrinsic colours to be (B–V)₀ = −0.075 ± 0.003 mag and (U–B)₀ = −0.228 ± 0.013 mag, with reddening E(B–V) = 0.026 ± 0.009 mag and extinction A_V = 0.085 ± 0.030 mag. Additionally, a fit of the Pecaut & Mamajek (2013) main sequence (MS) colour grid to the UBVJHKs photometry for HIP 81208, allowing the E(B–V) to vary, reveals that the reduced χ² < 1 fits result in estimates of T_eff = 10970 ± 350K, E(B–V) = 0.020 ± 0.013 mag and A_V=0.064± 0.041 mag for the target.

Both the reddening estimates derived above – using the Q-method and the spectral energy distribution (SED) fitting – rely heavily on the UBV photometry for a single star, namely HIP 81208, and are therefore not independent. However, both these estimates (E(B–V) ≤ 0.02±0.02 mag) are consistent with the STILISM reddening map from Lallement et al. (2018) (E(B–V) = 0.011±0.021 mag), which has the benefit of having already been published, and in addition is also a reddening estimate averaged over a number of stars of varying distances and positions in the vicinity of HIP 81208. The latter estimate is therefore more robust than the other estimates of reddening derived or cited in this section, although the consistency between them adds to its confidence. We therefore adopt the STILISM value of reddening and the inferred interstellar extinction of A_V = 0.034±0.064 mag for the target hereafter in this paper.

3.2 Spectral energy distribution and effective temperature

We used the Virtual Observatory SED analyzer (VOSA⁴) version 7.5 to query for the photometry of HIP 81208 and fit synthetic stellar spectra to the data. The observed spectrum for 78 photometric data points is shown in Fig. 1. We removed several discrepant or redundant photometry points, namely the synthetic Gaia DR3 data (for e.g. U, B, V, R, I) as they were dependent on other available Gaia DR3 photometry, and APASS Sloan photometry. Gaia DR3 photometry from binned B_p and R_p spectra account for 66 of the data points, along with the Gaia DR3 G_BP/G/G_RP points, and photometry from Tycho (BT, VT), APASS (B, V), 2MASS (J, H, K_s), WISE (W1, W2, W3, W4). For the SED fitting, we assumed a prior range of A_V = 0.034 ± 0.064 mag. Using Kurucz 2003 [a/Fe] = 0.0 models (Castelli & Kurucz 2003), the best overall fit that fits all 78 data points had the following parameters:

A_V = −0.03 mag, T_eff = 10750 ± 125 K, log(g) = 4, [^M/H] = −1, ℱ_bol = (8.43 ± 0.03) × 10⁻⁸ erg cm⁻² s⁻¹, and ℒ_bol = 58.27 ± 1.304 ℒ_⊙ (log(ℒ/ℒ_⊙) = 1.765 ± 0.01). The existing photometry as well as the obtained SED fit do not show any infrared excess up to ~10 µm (see lower panel of Fig. 1), which indicates the possible absence of any disk around HIP 81208. The SED fitting clearly favours models with negligible or zero reddening. The parameters all seem reasonable, except for the low metallicity, which was also found by the Gaia DR3 GSP Phot Aeneas fit to the Bp/Rp spectra ($\left[M/H\right]=-{1.27}_{-0.21}^{+0.39}$). Such a low metallicity is unexpected as the metallicity of Sco-Cen is essentially solar; however, the B9V star may be suffering from some abundance anomaly or other effect (e.g. fast rotation), which may be biasing the SED fitting.

We obtained a T_eff = 10 970 K from the SED fit of the MS colour grid to the UBVJHKs photometry in Sect. 3.1. From the inferred colours from the Q-method described in the same section, the latest version of the Pecaut & Mamajek (2013) table⁵ suggests that these values of intrinsic colours correspond to a T_eff of 11130± 110 K (and indeed very typical for B9V dwarfs).

Furthermore, the star’s spectral type itself (B9V) suggests a T_eff of approximately 10 700 K based on values for B9V standard stars by Pecaut & Mamajek (2013).

An additional independent T_eff estimate can be derived from the H_β photometry, β = 2.842±0.005 mag (see Slawson et al. 1992). Using the ATLAS9 model grid from Castelli & Kurucz (2003) for assumed atmospheric parameters log(g) = 4.0, [^M/H] =0, and υ_turb = 2 km s⁻¹, the Slawson H_β value corresponds to T_eff = 10 654 K (we note that adopting υ_turb = 0 km s⁻¹ has negligible effect, yielding T_eff = 10 679 K).

We therefore derived several T_eff estimates ranging from 10 654 K to 11 130 K. These five estimates are not quite fully independent, as three of the estimates are connected to the Pecaut & Mamajek (2013) dwarf colour/T_eff sequences, but they are consistent with an approximate T_eff of 10 840±220 K (rms). This value is within 2.5σ of the Gaia DR3 GSP-phot Aeneas estimate of $11{396}_{-76}^{+57}\text{K}$.

Fig. 1 Spectral energy distribution of HIP 81208 obtained with 78 existing photometry points retrieved from VOSA. Also shown is the best-fit stellar model (Kurucz ODFNEW /NOVER, alpha: 0.0 (2003); Castelli & Kurucz 2003) for a Teff = 10 750 K.

3.3 Luminosity and radius

We calculated a new bolometric luminosity for our target from the A_V and T_eff obtained in Sects. 3.1 and 3.2 and V mag from Table 1. We converted V to m_bol using a bolometric correction of BC_υ = −0.447 mag, obtained from interpolating the PM2013 table (refer footnote 5) for a T_eff = 10 840 K. The resulting apparent bolometric magnitude is m_bol = 6.154 ± 0.064mag. At a distance of 148.7 pc, this corresponds to the absolute magnitude M_bol = 0.292 ± 0.080 mag. This translates to an equivalent luminosity of ℒ_bol = 60.69 ± 4.46 ℒ_⊙ (log(ℒ/ℒ_⊙) = 1.783 ± 0.032), which is consistent with previous ℒ_bol values for HIP 81208 in the literature (e.g. Anderson & Francis 2012).

Using the above derived value for ℒ_bol and the T_eff estimate from Sect. 3.2 (10 840 ± 220 K), we can arrive at an estimate for the radius of the target using the Stefan-Boltzmann law, which gives ${ℛ}_{*}=\sqrt{{ℒ}_{\text{bol}}/4\pi \sigma T_{\text{eff}}^4}=2.213\pm 0.121R_{\odot }$.

3.4 Stellar age

The stellar age was constrained using the method of isochronal dating of co-moving stars (CMS) introduced in the previous BEAST publications, namely Janson et al. (2021a) and Squicciarini et al. (2022). We adopted the long-term proper motion for HIP 81208 reconstructed by Kervella et al. (2022) for this purpose. The sample of co-moving stars was constructed by querying the Gaia DR3 catalogue within a search radius of 5 deg around the location of our target, for sources possessing similar 2D-projected space motion. A clump of 86 such sources was identified in the velocity space (υ_α, υ_δ) through the following cuts on velocity and parallax (ϖ): $$\{\begin{array}{l}-7.5\text{\hspace{0.17em}\hspace{0.17em}km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}<{\upsilon }_{\alpha }<-5.8\text{\hspace{0.17em}\hspace{0.17em}km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}\hfill \\ -19.5\text{\hspace{0.17em}\hspace{0.17em}km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}<{\upsilon }_{\delta }<-15.5\text{\hspace{0.17em}\hspace{0.17em}km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}\hfill \\ 6.0\text{\hspace{0.17em}\hspace{0.17em}mas\hspace{0.17em}\hspace{0.17em}<\hspace{0.17em}}\varpi <7.3\text{\hspace{0.17em}\hspace{0.17em}mas.}\hfill \end{array}$$(1)

Based on their photometry from Gaia DR3 (G, G_BP, G_RP) and 2MASS (J, H, K_s) bands, we subsequently derived the isochronal ages and masses for all the stars in the group using the MADYS tool (Squicciarini & Bonavita 2022). Three independent evolutionary models, all assuming solar metallicities, were employed to ensure the robustness of the results: BHAC15 (Baraffe et al. 2015), PARSEC (Marigo et al. 2017), and MIST (Choi et al. 2016). The sample was then divided into five mass bins: ^M/M_⊙≤0.4, 0.4<M/M_⊙≤0.6, 0.6<M/M_⊙≤0.8, 0.8<M/M_⊙≤1.4, ^M/M_⊙>1.4. We restricted our age analysis to the most reliable bin of stellar mass (0.8 M_⊙ < M ≤ 1.4 M_⊙), which contains 9 out of the 86 stars in the sample. The small number of stars in this bin is a consequence of the shape of the initial mass function (IMF; 53 of the fitted stars were M-type stars, i.e. M < 0.6 M_⊙) and the fact that not all the stars could be successfully fitted because of factors such as non-optimal quality of photometric data or unresolved binarity. Using the chosen mass bin, we were able to derive a well-defined group age estimate $t={17}_{-4}^{+3}\text{Myr}$, which is consistent with previous age estimates for UCL, such as those of Mamajek et al. (2002, 17±1 Myr) and Pecaut et al. (2012, 16±1 Myr).

3.5 Stellar mass and log(g)

To estimate the stellar mass, we used the Parsec tracks (Bressan et al. 2012) version 1.2S (Chen et al. 2015) with A_V = 0.034 mag and metallicity [^M/H] = 0 and age = 17 Myr. We interpolated for T_eff = 10840±220 K and log(ℒ/ℒ_⊙) = 1.7831 ±0.0319 to get an estimate for the mass of the star, ℳ_* = 2.58±0.06 ℳ_⊙, and its surface gravity, log(g) = 4.201±0.011. Uncertainties for ℳ_* and log(ɡ) were derived using the upper and lower limits for T_eff(±220 K) and log(ℒ/ℒ_⊙) (±0.0319). To get the upper limit on the mass and surface gravity, we interpolated the Parsec tracks for ℳ_* and log(g) corresponding to T_eff = 11060 K and log(ℒ/ℒ_⊙) = 1.815. Similarly, to get the lower limits, we interpolated using T_eff = 10620 K and log(ℒ/ℒ_⊙) = 1.7512. The average of the uncertainties derived from the respective upper and lower limits of ℳ_* and log(g) obtained in this way are quoted as the uncertainty for these parameters in this section. The values for ℳ_* and log(ɡ) derived here are also consistent with the literature (see e.g. Anders et al. 2022).

4 Results and discussion

Post-processing of the IFS image using ADI reveals a candidate companion, HIP 81208B, at a projected separation of 0.325 ± 0.001′′ (48.33 ± 0.48 au; averaged over the two epochs) north of the primary star. The candidate is almost equally bright in both the epochs (ΔJ = 7.6 ± 0.1 mag), as can be seen in Figs. 2a and 2b, which display the median of the reduced images over the entire IFS YJH-band in the 2019 and 2022 epochs, respectively. The candidate is also detected in the reduced data from IRDIS K1 (ΔK1 = 6.9 ± 0.1 mag) and K2 (ΔK2 = 6.6 ± 0.1 mag) bands in both epochs (refer Figs. 3a and b). In addition, a second candidate companion, HIP 81208C, is also revealed in the IRDIS images at a projected separation of 1.492 ± 0.001″ (221.86 ± 2.09 au, averaged over the two epochs) from the central star. This candidate is not within the field of view (FoV) of the IFS (1.73″ × 1.73″).

Table 3 lists the obtained photometry and astrometry for HIP 81208 B and C from the observations in both epochs. To put both candidates in the same astrometric reference frame, the astrometry was calculated based only on the IRDIS data. The astrometry in the individual IRDIS bands is obtained by inserting a negative point spread function (PSF) in the approximate position of the planet, and minimising the root mean square (rms) of the residuals in an area of ~1λ/D around this position in the reduced data by appropriate changes in both the exact position where the PSF was inserted and the multiplicative factor needed to reproduce its intensity. The final astrometry values listed for the companion are the weighted mean of its values from K1 and K2 bands. The contrast magnitudes for the companions were obtained using a similar procedure, but keeping the position of the companion fixed at the value obtained from the astrometric analysis. For IFS photometry, all the individual 39 channels were considered separately (see Sect. 4.3 for more details) and the contrast magnitudes in Y, J, and H bands were obtained by taking the weighted mean of the respective values over the spectral channels with wavelength between 1–1.1 µm, 1.16–1.33 µm, and 1.5–1.64 µm, respectively. Figure 4 shows the position of the two candidate companions in the K1 – K2, K1 colour–magnitude diagram, along with other known substellar companions, young and/or dusty objects, and some field stars.

4.1 Confirming the physically bound nature of the candidates

In order to confirm whether or not the directly imaged candidates are physical companions to the primary, we study their proper motion (see e.g. Janson et al. 2021a). A physical companion will have a similar proper motion to the star, but one that is different from those of other background sources in the FoV. Aside from the primary star and the two companion candidates HIP 81208 B and C, we identify six other sources detected within the IRDIS (K1, K2) FoV that are common to the data in both epochs. Five additional sources are also seen, but only in one epoch, and these are therefore neglected in the following analysis. We study the astrometric motion of all of the above objects relative to the star between the two epochs (given in Table A.1, along with the obtained photometry), as shown in Fig. 5. Given the magnitude of the two candidates B and C, if they are not physically bound to the target star, they would need to be background stars and would in that case possess a low proper motion tending towards zero, that is, a null proper motion. Their astro-metric positions would then shift between the two epochs as a reflection of the proper motion of the star. This astrometric motion can be predicted and is shown as the black dotted curve in the figure, starting at the position expected for a source co-moving with HIP 81208 (unfilled star) and ending at the position expected from the astrometric shift of a pure background source with null proper motion in the second epoch (black filled star). Also plotted in the figure is the relative astrometric shift of the six additional sources detected in the IRDIS data in both epochs (blue crosses, labelled as ‘prob. bkg objects’) along with the associated error bars. These shifts are large with respect to the predicted motion of a comoving source; the offsets from the position expected for a completely static source are due to their own proper motion, which could be the signature of a reflection of the galactic rotation curve. It can be seen that the motions of the two companion candidates B (green cross) and C (yellow cross) are clearly different from those of these six sources, and are much closer to that of a comoving source. The mean astrometric shift of the six additional sources in the IRDIS data is −12.7 ± 1.5 mas (rms: 5.1 mas) along αcos(δ) and −14.7 ± 1.3 mas (rms: 4.4 mas) along δ, plotted as the thicker blue cross; on the other hand the astrometric shifts along αcos(δ) and δ of HIP 81208 B and C are (7.6 ± 9.6, −62.7 ± 1.4) mas and (−12.9 ± 2.2, −62.5 ± 2.2) mas, respectively. This correspondingly places a confident 10.4 σ and 9.8 σ distinction between the astrometric motion of B and C from the mean astrometric motion of the cloud of background objects detected in the IRDIS field, strongly indicating that these candidates are physically bound to the primary star.

To quantify this indication, we compute the probability that the two companion candidates, if drawn from such a sample of background sources, could possess an astrometric shift that would mimic that of a physically bound companion of the primary. In order to compute such a false alarm probability (FAP), we queried the Gaia EDR3 catalogue (Gaia Collaboration 2020) for field stars near HIP 81208 with a search radius of 5 arcmin. The resulting sample contains 4014 Gaia sources; the 50th percentile of their proper motion distribution (with the 16th and 84th percentiles shown as the lower and upper uncertainty limits) is $-{3.12}_{-2.71}^{+3.11}$ mas yr⁻¹ along α^* (〈µ_α*〉) and $-{3.67}_{-2.85}^{+2.57}$ mas yr⁻¹ along δ (〈µ_δ〉). These sources are shown as the grey dots in Fig. 6, labelled ‘Gaia sources’. In the same figure, we also show the six background sources from IRDIS data as blue crosses, with their proper motion computed as the ratio between the corresponding astrometric shifts with respect to the position expected for a static source and the time baseline between the two epochs. The mean proper motion of these six sources, along with the associated uncertainty, is 〈µ_α*〉 = −4.77 ± 0.58 mas yr⁻¹ and 〈µ_δ〉 = −5.50 ± 0.50 mas yr⁻¹ (shown as the thicker blue cross in the figure), which is well within ~1σ of that of the Gaia sample. Similarly, the derived µ_α* and µ_δ of HIP 81208 B and C are (2.87 ± 3.62 and −23.53 ± 0.54) mas yr⁻¹and (−4.85 ± 0.84 and −23.46 ± 0.75) mas yr⁻¹, respectively. To estimate the FAP, we define an ‘interesting region’ of proper motion in the figure, similar to Squicciarini et al. (2022), within which a background star might be disguised as a physically bound source. This region employs the following boundaries: $$\{\begin{array}{l}-25\text{\hspace{0.17em}\hspace{0.17em}mas\hspace{0.17em}\hspace{0.17em}}{\text{yr}}^{-1}<{\mu }_{{\alpha }^{*}}<5\text{\hspace{0.17em}\hspace{0.17em}mas\hspace{0.17em}\hspace{0.17em}}{\text{yr}}^{-1},\hfill \\ {\mu }_{\delta }<-20\text{\hspace{0.17em}\hspace{0.17em}mas\hspace{0.17em}\hspace{0.17em}}{\text{yr}}^{-1}.\hfill \end{array}$$(2)

The number of Gaia stars that fall into this region is 12 out of 4014 stars, resulting in a corresponding fraction of ‘interesting objects’ of ~3 × 10⁻³. Given that there are eight source detections in the entire IRDIS FoV, this gives a FAP of ~2%. The probability of having at least two such sources is as low as 2 × 10⁻⁴. The obtained low value of FAP is sufficient to add further robustness to the claim that HIP 81208 B and C are physically bound companions to the primary.

It is also evident from Table A.1 that both the companions B and C are much brighter than all other sources detected in the IRDIS FoV. An average over the K1 and K2 bands over the two epochs would yield K = 13.51 ± 0.06 mag for B and K = 12. 44 ± 0. 06 mag for C. An alternative way to determine if HIP 81208 B and C are background objects would then be to compute the probability that background objects as bright as these can be detected at such small separations from the star. As the two companions are quite bright, we may use the 2MASS catalogue (Skrutskie et al. 2006) – which is complete for uncrowded regions at this magnitude level – to estimate the surface density of background objects as bright as B and C. Searching within a radius of 2 arcmin from HIP 81208, we found 26 and 15 2MASS sources brighter than B and C, respectively. Given the projected separation of the two companions from the star, this corresponds to very low probabilities of respectively 1.9 × 10⁻⁴ and 2.3 × 10⁻³ that B and C are background objects by chance projected very close to HIP 81208.

An underlying possibility that is yet to be investigated is that the two companions, although co-moving as established above, could be UCL members that are not bound to the target. Considering the derived masses of B and C, we first evaluate the probability that B is a free-floating UCL substellar object that was by chance close to the target. For this, we perform a similar analysis to that of Squicciarini et al. (2022) and integrate the normalised IMF of the UCL association from 5 M_J⁶ to 75 M_J to obtain the fraction of objects in UCL belonging to this mass range. To obtain the expected number of UCL objects in the range of 5–75 M_J, we multiply this fraction by the number of UCL sources (~4021)⁷. We further divide this value by the area of the UCL subregion as given in Squicciarini et al. (2022) to estimate the corresponding projected density of these objects. Multiplying the projected density with the IRDIS FoV (11″×11″) gives us the number of free-floating substellar objects expected in the IRDIS FoV. Using a binomial distribution, we obtain a FAP of having seen at least one such object across the 47 targets that have been observed at least twice in the BEAST survey of ~5 × 10⁻⁴. Similarly, we also computed the probability that C is a stellar mass UCL member that is not bound to HIP 81208 by repeating the above calculation but considering every possible companion mass from 5 M_J to the most massive UCL star ~15 M_⊙. The corresponding value for FAP is ~2.4 × 10⁻³. This probability can be regarded as an upper limit, because this analysis considers all the stars in UCL, including those for which Gaia has resolved the true position. For a more accurate estimation, we should consider only those stars within the IRDIS FoV that are undetected by Gaia (e.g. because of being hidden behind brighter stars); this would mean the sample of interloping stars from which the FAP is computed would be smaller than the current one. Therefore, we can safely say that the probability of seeing at least one UCL member of any stellar mass >5 M_J across all the stars that have been observed twice in the BEAST survey until now by chance due to projection effects is <0.2–0.3%.

Fig. 2 SPHERE IFS YJH band image of HIP 81208 from (a) the 2019 epoch and (b) the 2022 epoch, reduced using a KLIP-based pipeline. The image here is the median taken over the entire wavelength band and shows a bright companion, HIP 81208 B, detected at 0.325 ± 0.001″ north of the primary star (averaged over the two epochs). The pixel scale is 7.46 ± 0.02 mas pixel−1 in both epochs.

Fig. 3 SPHERE IRDIS K band image of HIP 81208 from the (a) 2019 and (b) 2022 epoch observations, reduced using SpeCal in TLOCI mode. Each image here is the median taken over K1 and K2 band images of the respective epoch. Images show the companion HIP 81208 B, as well as a second companion HIP 81208 C detected at 0.325 ± 0.001″ north and 1.492 ± 0.001″ southeast of the primary star, respectively (averaged over the two epochs). The pixel scale was 12.258 ± 0.004 mas pixel−1 for K1 band and 12.253 ± 0.003 mas pixel−1 for K2 band.

Fig. 4 K1 – K2, K1 colour–magnitude diagram, with the candidates HIP 81208B and HIP 81208C shown along with other known substellar companions, young and/or dusty objects, and field stars.

Fig. 5 Relative astrometric shifts between the two epochs along α* (x-axis) and δ (y-axis) for all the sources detected in the IRDIS FoV. The black dotted curve shows the expected motion of a pure background source with null proper motion when compared to a source co-moving with the target (the grey shaded region encompasses the uncertainty associated with this motion). The mean astrometric shifts of the probable background objects (blue crosses) detected in the data are shown as a thick blue cross, and those of HIP 81208 B and C are shown as a solid green triangle and a solid yellow square, respectively.

Fig. 6 Proper motion in αcos(δ) and δ for the sample of Gaia EDR3 sources within 5 arcmin of the primary, and for the six background sources in the IRDIS data. Also shown are the median proper motion of the six probable background objects as well as the proper motion derived for HIP 81208 B and C using the same symbols as in Fig. 5. The interesting region of proper motion defined by Eq. (2) is over-plotted in the figure as a yellow-shaded region enclosed within black dashed contours.

4.2 Characterisation of candidate properties

Photometric mass estimates were derived for both the candidates with the help of MADYS based on the average of their respective contrast measurements over the two epochs (shown in Table 3). In particular, the estimate for HIP 81208C is based on IRDIS (K1, K2) bands while that for HIP 81208B is based on both IRDIS and synthetic IFS magnitudes J_IFS at 1.246 µm (band width = 0.174 µm) and H_IFS at 1.570 µm (band width = 0.132 µm), which were obtained by collapsing the spectral channels 12–21 (1.159–1.333 µm) and 30–38 (1.504–1.636 µm), respectively. The derivation of calibrated apparent magnitudes from contrasts is mediated by suitable photometric data for the primary taken from the literature: 2MASS J, H, and K_S for IFS J_IFS, H_IFS, and the doublet (K1, K2), respectively (shown in Table 1). Given the spectral classification of the primary as a B9V star, the errors introduced by these approximations are well within photometric uncertainties.

MADYS operates by seeking the best match between the vector of input photometry for the object and a selected grid of isochrone and evolutionary track via χ² minimisation. Evaluation of random uncertainties associated to apparent photometry, parallax, extinction, and age is naturally taken into account in this process in a Monte Carlo fashion. In order to evaluate the impact of theoretical uncertainties on the results, we repeated the mass estimation of each companion twice, employing a different suite of stellar and substellar evolutionary models each time; namely the Ames-Cond models (Allard et al. 2001) and the BHAC15 models. For each candidate, the individual estimates from both repetitions were consistent with one another, and were averaged to get the final mass estimate: $$M_{\text{B}}={67}_{-7}^{+6}{M}_{\text{J}}={0.064}_{-0.007}^{+0.006}{M}_{\odot },$$(3) $$M_{\text{C}}={141}_{-14}^{+10}{M}_{\text{J}}={0.135}_{-0.013}^{+0.010}{M}_{\odot }.$$(4)

The corresponding T_eff returned by the method as a result of the above interpolation scheme were also averaged to get $\langle T_{\text{eff}}\rangle ={2895}_{-40}^{+45}\text{K}$ for B and $\langle T_{\text{eff}}\rangle ={3165}_{-60}^{+40}\text{K}$ for C. As indicated from Table 6 of Pecaut & Mamajek (2013), which lists intrinsic colours, adopted T_eff, and bolometric corrections of 5–30 Myr stars, these values roughly correspond to a spectral type of M5 for B and M4 for C. A comparison of the T_eff estimate for B can be made using its observed spectrum over the IFS and IRDIS bands (see Sect. 4.3), giving us some level of confidence as to the derived T_eff and spectral type. However, no such comparison can be made for C because it lacks an observed spectrum in the IFS, and so the derived T_eff and spectral type for C here should be subjected to careful interpretation.

Grieves et al. (2021) briefly summarises the stellar–substellar mass boundary predictions in the literature. Although the hydrogen-burning mass limit is generally adopted as 80 M_J, the boundary predictions range from 73.3 to 98.5 M_J depending on the choice of models and metallicities (see Dieterich et al. 2018). An even lower mass limit of ~70 M_J was predicted by Dupuy & Liu (2017) based on the astrometric masses of 31 ultracool binaries with component spectral types M7–T5. Comparing the above derived mass estimates for the companions with the stellar–substellar boundary predictions will place HIP 81208C in the stellar regime and HIP 81208B most likely in the brown dwarf regime. However, further observations and more stringent mass constraints are required to definitively confirm whether HIP 81208B is a brown dwarf or a very low mass star.

4.3 Analysis of HIP 81208B spectrum

To determine contrast values for HIP 81208B from the KLIP-reduced IFS data, a negative PSF was injected at the position of the candidate. Monochromatic reductions were then performed based on principal component analysis (PCA) with modes 2–6 in such a way that the standard deviation is minimised within a region of ~1λ/D centred on the position of the candidate. The values obtained for each wavelength with the different PCA modes were averaged to get the final contrast for that wavelength. The uncertainty in the contrast for each wavelength was calculated as the rms scatter of the values obtained by repeating the above same procedure, with the negative PSF inserted, not at the companion position but at five different positions in the image at the same separation from the star as the companion B, and separated by 60, 120, 240, and 300 degrees in position angle. The minimisation of the standard deviation only concerns the intensity of the negative PSF that is inserted and can therefore be performed even without an actual companion present at the location of insertion. The intensity value derived this way can be either positive or negative, depending on the local noise realisation. The contrast values and uncertainties thus obtained for the candidate from the IFS data in both the epochs are plotted against the respective wavelengths in Fig. 7 along with the obtained contrasts using SpeCal from the IRDIS K1 and K2 bands.

To calculate the flux of the companion B corresponding to the contrast values at each wavelength, we interpolated the bestfit theoretical spectrum for the star from Sect. 3.2 (see Fig. 1) scaled to the observed flux levels at the respective wavelengths. The uncertainties in flux are obtained by quadratically summing the uncertainties in contrast and the photon noise at the location of B estimated within an area of ~1λ/D around it. Further, to determine the spectral type of the candidate, we fit theoretical stellar models AMES-Cond 2000 (Allard et al. 2001) and BT-Settl (Allard et al. 2012) to the obtained flux spectrum of the candidate. To determine the right choice for surface gravity, we interpolated the AMES-Cond isochrone corresponding to the estimated age of the target in this work (0.02 Gyr) for the estimated mass of HIP 81208 B in Sect. 4.2 (0.064 M_⊙); the corresponding log(g) given by the isochrone is 4.42. The closest available choice of surface gravity among the models is log(ɡ) = 4.5. As an additional check, we also interpolated the AMES-Cond evolutionary track corresponding to a mass of 0.06 M_⊙; at a log(ɡ) = 4.5 the age suggested by the evolutionary track is 0.026 Gyr, which is compatible with our age estimate for the target. For these reasons, we chose a log(ɡ) value of 4.5 for the theoretical models in this analysis. The resulting theoretical flux at the surface of the star was scaled to the observed flux level by a factor of R²/D², where R is the radius of the candidate for the specific temperature of the model as determined from the corresponding model isochrones and evolutionary tracks corresponding to an age of 0.02 Gyr, and D is the distance to HIP 81208 from Earth. Figure 8 shows the reduced χ² (${\chi }_{\text{red}}^2$) values between the observed spectrum and the two theoretical models at specific effective temperatures. For both 2019 and 2022 epoch observations, the ${\chi }_{\text{red}}^2$ value for both models approaches the ideal value of 1 at the effective temperature of 2900 K, the closest spectral type to which is M5 (T_eff = 2880 K) (see Table 6 of Pecaut & Mamajek 2013). Figure 9 shows the observed spectrum for HIP 81208B in both the epochs alongside theoretical models corresponding to the lowest ${\chi }_{\text{red}}^2$ (2900 K), as well as for T_eff = (2800 K, 3000 K). The obtained T_eff here for B from its spectrum and the corresponding spectral type are consistent with the results from MADYS in Sect. 4.2, supporting the reliability of these values. We therefore conclude that HIP 81208B is likely of spectral type M5.

Fig. 7 Contrast spectrum for HIP 81208B obtained from reduced IFS data in the YJH band and IRDIS data in K1 and K2 bands, along with the associated uncertainties. The red solid squares represent the contrast obtained from the 2022 observations and the green solid circles represent the contrast obtained from 2019 observations.

Fig. 8 ${\chi }_{\text{red}}^2$ computed between the observed spectrum and the AMES-Cond and BT-Settl models for different effective temperatures for both 2019 and 2022 epochs. The log(ɡ) adopted for both models is 4.5. For both epochs, the models reach a ${\chi }_{\text{red}}^2\sim 1$ at Teff = 2900 K. The closest spectral type is M5.

4.4 Constraining the orbital parameters

Given the astrometric information of the two companions obtained in two different epochs, it is possible to derive constraints on their orbital parameters. We do this using the Python package orbitize! (Blunt et al. 2020), which is designed to fit the orbits of directly imaged planets. Orbitize! offers two choices of algorithm for fitting orbits: Orbits for the Impatient (OFTI, see Blunt et al. 2017) and Markov chain Monte Carlo (MCMC, see Ford 2005; Vousden et al. 2016). The disadvantage with using MCMC for orbital fitting of long-period orbits is that for cases with less constraints from observations on the orbital parameters, MCMC take a very long time to converge. As the relative displacement between the two epochs for B and C were only 7.8 mas and 3.4 mas, respectively, we had very low orbital coverage for the companions from our observations, and therefore MCMC would not be an ideal choice. This means that our choice of algorithm was OFTI, which is a Bayesian Monte Carlo rejection-sampling method that is ideal for cases where the observations cover only a small fraction of long-period orbits. OFTI takes the separation and position angle of the companion in the different observation epochs, and the parallax and the total mass of the system (star and companion), along with uncertainties, as input and uses built-in⁸ prior probability distribution functions (PDFs) to compute posterior PDFs of orbital parameters. This computation is based on the orbits that the algorithm accepts from the generated ones, using the technique of rejection-sampling. We ran orbitize! until a total of 10⁶ orbits were accepted by the algorithm for both B and C. The corner plots in Figs. 10a and b show the resulting posterior distribution of the semi-major axis (a), inclination (i), and eccentricity (e) of these orbits (the plots are truncated to include only 98% of the generated orbits for better visibility). The 2D contour plots in the figures illustrate 1σ, 2σ, and 3σ confidence levels on these values. The median orbital parameters for B and C as derived from these posteriors, along with the upper (84%) and lower (16%) quantile intervals (~±1σ), are also shown in the figures, and are listed in Table 4. The orbital periods implied for B and C from their semi-major axis obtained above and their derived masses in Sect. 4.2 are also listed in the table.

The corner plots suggest it is very unlikely that B has a very high eccentricity orbit (>0.7); higher inclinations (>75º) than the predicted quantile limit would favour a much wider orbit, both of which seem of very low probability based on the posteriors. Similarly, for C as well, the probability falls rapidly for higher eccentricities and higher values of orbital distance seem very unlikely. It is noteworthy from these results that the orbital inclination of C is very different from that of B, with the relative inclination between the two companions being ~80º (refer Fig. 11 for the distribution of relative inclinations between the orbits of B and C obtained from the algorithm). Thus, the two orbital planes look roughly orthogonal, making B and C appear to orbit in opposite directions.

Figures 12 and 13 show 20 random orbits drawn from the posterior distributions from the OFTI run for B and C, respectively. The left panels show the orbital motion of the respective companions in RA and DEC over an entire period of these orbits and the right panels show the projected separation ρ (mas) and position angle (PA; in deg) predicted from these orbits over time. Also shown in the plots are the observed ρ and PA of the companions in the two epochs along with error bars.

Fig. 9 Observed spectrum for HIP 81208B in 2019 (green solid circles) and 2022 (red solid squares) epochs, along with the associated uncertainties in flux. Also plotted are the AMES-Cond 2000 log(g) = 4.5 theoretical model for Teff = 2800, 2900, and 3000 K and the BT-Settl log(ɡ) = 4.5 theoretical model for Teff = 2900 K. The best-fit models to the observed spectrum in both epochs correspond to Teff = 2900 K.

Fig. 10 Corner plots for the accepted orbits from OFTI for companions B (a) and C (b). The median values along with upper (84%) and lower (16%) quantile intervals (~±1σ) are shown as dashed vertical lines on the posterior distributions of the respective orbital parameters. The 2D contour plots show the 1σ, 2σ, and 3σ confidence limits to these values.

4.5 Dynamical stability of the orbits

We can refine the current orbital constraints (depicted in Fig. 10) by ensuring that the configuration is stable for a rimescale comparable to the age of the system (~10⁷ yr). Indeed, it is unlikely that we observe this system immediately before the ejection of one of its components. The stability of a triple system with non-negligible mass ratios requires a well-separated hierarchical structure, with sufficient distance between the orbits. Moreover, if the eccentricities and/or relative inclination are high, then additional secular perturbations can destabilise orbits that were initially well separated. To evaluate the additional constraints ensuing from dynamical stability, we randomly pick 1000 couples of solutions to the orbital fits, and run N-body simulations of their future 10⁷ yr of evolution. We use the IAS 15 integrator in the Rebound package (Rein & Spiegel 2015). We adopt the following values for the masses: M_* = 2.58 M_⊙, M_B = 0.064 M_⊙, and M_C = 0.135 M_⊙. The maximum value of the eccentricity reached by either orbit is displayed in Fig. 14 as a function of the semi-major axes, eccentricities, and relative inclination. Values of more than 1 signify instability (corresponding orbits are indicated as black solid circles in the figure). The results suggest that a large fraction of the orbital solutions (nearly half of them) are actually unstable. In particular, for companion C, the lower right panel of Fig. 14 shows that instability develops starting from e_c ≃ 0.5, whatever the relative inclination (the chaos becomes too strong for the system to be locked in a Kozai regime). From the top-middle panel of the figure, we also see that only wide separations (a_c > 500 au) are able to help maintain the stability for e_c > 0.5; however, such separations are outside the interval of confidence for C (see Table 4). Hence, we can say that large eccentricities (>0.5) are disfavoured for C.

In the simulations, we notice that a significant proportion of solutions reach high eccentricities without becoming unstable. This is characteristic of the Von Zeipel-Lidov-Kozai (ZLK Von Ziepel 1910; Lidov 1962; Kozai 1962) resonance, which occurs when the relative inclination is more than ≃40 degrees and less than ≃130 degrees. Looking at the whole set of solutions, the likelihood that this system is in a configuration of ZLK resonance is around 65% (see Fig. 11). The period of Kozai eccentricity oscillations for a triple system can be approximately estimated as (Takeda & Rasio 2005; Ford et al. 2000): $$P_{\text{Koz}}\simeq P_1\left(\frac{m_0+m_1}{m_2}\right){\left(\frac{a_2}{a_3}\right)}^3{\left(1-e_2^2\right)}^{3/2},$$(5)

where P, a, and e are the orbital period, semi-major axis, and eccentricity, respectively, and the indices 0, 1, and 2 are in this case representative of the primary star HIP 81208, the brown dwarf companion B, and the stellar companion C, respectively. For the system HIP 81208, using the derived mass and orbital parameters in Sects. 4.2 and 4.4, the Kozai oscillation period is ≃0.3 Myr, which is much shorter than 17 Myr, which is the age of the system. Therefore, the system has enough time to have undergone at least a few cycles of Kozai oscillations. Such systems are rare and precious, for ZLK resonances are at the core of many astrophysical problems such as hot Jupiters (see Wu & Murray 2003) or black hole mergers (see Liu & Lai 2019). More generally, the relative inclination between the different orbits of a given system is far from well constrained, and could bring important clues as to the formation processes. In this regard, HIP 81208 probes an uncharted regime between the planetary (Dupuy et al. 2022) and stellar masses (Tokovinin 2017), where it is unclear whether coplanarity or random inclinations should be favoured.

Fig. 11 Distribution of the relative inclinations between B and C’s orbits from the sample of solutions obtained by the OFTI algorithm. The light red area corresponds to systems in the ZLK configuration, or Kozai regime.

4.6 A potential ‘D’ component

In addition to the two companions reported in this study, the star Gaia DR3 6020420074469092608 (2MASS J16360769-3543514, WISEA J163607.70-354351) at a separation of 656.28″ (0.1823°) and position angle of 92.046° is potentially noteworthy. The star is a faint (G = 14.975 ± 0.003 mag; Gaia Collaboration 2020) X-ray-emitting source detected by the Rosat All-Sky Survey (leRASS J163607.7-354351, 2RXS J163607.8-354352; see Boller et al. 2016). Its parallax, 6.6742 ± 0.0297 mas (Gaia Collaboration 2020), differs from that of HIP 81208 only by ~3σ and the proper motion (µ_α* = −9.622 ± 0.035 mas yr⁻¹, µ_δ = −24.695 ± 0.029 mas yr⁻¹; Gaia Collaboration 2020) only by ~ 1.12 mas yr⁻¹(856 ± 56 ms⁻¹). Given the total mass of the HIP 81208 triple system, namely M_tot = 2.779 M_⊙, the tidal radius of HIP 81208 is then 1.35 pc×(M_tot/M_⊙)^1/3 = 1.898 pc (Mamajek et al. 2013; Jiang & Tremaine 2010). At a distance of 148.7 pc, this translates to a projected tidal radius of (180/π)×(1.898 pc/148.7 pc) = 0.7313°, or 2632.75″. As the star is within 656.28″of HIP 81208, it is possible that it could be within the tidal influence of our target. This fact, along with the similar parallax and proper motion to HIP 81208, suggest a possibility that this star could be a potential ‘D’ companion in this currently triple system. However, we cannot reliably assess this possibility given the present lack of knowledge about this star. Future observations of Gaia DR3 6020420074469092608 with instruments such as GRAVITY (GRAVITY Collaboration 2017) could help provide more information on its projected 2D motion in the sky. In addition, radial velocity information from future spectroscopic observations will also help determine its 3D motion. If the star shares the same 3D motion in the sky as HIP 81208, this could be a strong indication that it is bound to our target.

5 Conclusions

In this work, we describe our observations of the B9V star HIP 81208 using the SPHERE instrument at VLT as part of the BEAST survey in two different epochs, obtaining both IFS YJH spectroscopy and IRDIS K1 and K2 imaging data for the target. We report the discovery of two lower-mass companions to the star; an inner companion HIP 81208B at a projected separation of 0.325 ± 0.001″ and an outer companion HIP 81208C at a projected separation of 1.492 ± 0.001″. Using its spectrum over the YJH – K1 – K2 wavelength range, we determine that HIP 81208B is of spectral type M5. Analysing the obtained photometry, we estimate the masses of HIP 81208B and C to be ${67}_{-7}^{+6}{M}_J$ and ${0.135}_{-0.013}^{+0.010}{M}_{\odot }$, respectively, indicating that B is most likely a brown dwarf and C a low-mass star. We confirm that both companions are physically bound to the primary by means of proper motion analysis of the background stars in the data.

Using astrometric information over the two epochs, we were able to constrain the orbits of B and C to the most probable semi-major axis values of ${53.98}_{-15.00}^{+32.22}$ au and ${234.27}_{-68.96}^{+168.65}$ au, respectively. The relative inclination between the two orbits is high, making them appear to orbit in an opposite sense to each other. Table 4 summarises the main parameters we derived for the two companions in this work. The orbital solutions we derived for the companions indicate that this system is likely to be in a Kozai resonance, in which orbits with high relative inclinations can remain stable despite reaching high values of eccentricity. This means that HIP 81208 is a system of considerable interest; careful further study of this system would help us to understand more about orbital dynamics in systems where the companions are in between planetary and stellar mass regimes.

Our observations therefore provide strong evidence that HIP 81208 is a triple system, with two low-mass companions in a very interesting configuration around the star. If the nearby star Gaia DR3 6020420074469092608 turns out to be physically bound to HIP 81208, the system may even be quadruple. Further observations of HIP 81208 would therefore be of use, both in investigating the possibility of perhaps a ‘D’ component in the now known triple system, and in understanding more about the orbital dynamics in the system. The orbital parameters derived in this study for HIP 81208 B and C are not tightly constrained – as can be seen from their large ±1σ errors – because they are based on a small time baseline. Future BEAST or GRAVITY epochs will aid in providing longer time baseline astrometry for both the companions, which could help establish more stringent constraints on the orbital parameters for this very interesting system.

Fig. 12 Illustration of 20 random orbits drawn from the posterior distributions from the OFTI run for HIP 81208B. The left panel shows the motion of B in RA and Dec over an entire orbital period for these 20 cases. The adjacent color bar shows the epochs along the orbital positions. The yellow solid circles represent the predicted position of the companion in these orbits in the two observation epochs. The blue star represents the central primary star HIP 81208. The right panel shows the projected separation ρ and PA (deg) of B over time, as predicted from these 20 orbits. The two purple dots on both figures represent the measured values of the respective parameters at the two observation epochs.

Fig. 13 Illustration of 20 random orbits drawn from the posterior distributions from the OFTI run for HIP 81208C. The left panel shows the motion of C in RA and Dec over an entire orbital period for these 20 cases. The change in the position of the companion over the two epochs is very small and the separate data points are not visible as in the previous figure and are therefore represented by a single yellow solid circle. The right panel shows the projected separation ρ and PA (deg) of C over time as predicted from these 20 orbits. The symbols and colour bar hold the same meaning as in the previous figure.

Fig. 14 Maximum eccentricity reached by either orbit within 107 yr of evolution depending on the initial orbital elements. Black solid circles depict unstable cases (emax > 1, 44% of the cases). These results were obtained by running N-body simulations starting from 1000 solutions of the orbital fit. Configurations where orbits are crossing or almost crossing are always unstable. Orbits with relative inclination between 40 and 130 degrees are locked in a Kozai resonance, and may experience high eccentricity variations without compromising the stability of the system.

Acknowledgements

M.J. gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation. A.V. and P.D. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 757561, COBREX; grant agreement no. 885593). This publication makes use of CDS and NASA/ADS services, as well as VOSA, developed under the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/through grant PID2020-112949GB-I00. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement no. 776403 (EXOPLANETS-A). This work has also 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). This research has made use of the SIMBAD database and VizieR catalogue access tool, operated at 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 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 is also supported by the PRIN-INAF 2019 “Planetary systems at young ages (PLATEA)”. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

Appendix A Properties of IRDIS sources

References

All Tables

All Figures

	Fig. 1 Spectral energy distribution of HIP 81208 obtained with 78 existing photometry points retrieved from VOSA. Also shown is the best-fit stellar model (Kurucz ODFNEW /NOVER, alpha: 0.0 (2003); Castelli & Kurucz 2003) for a Teff = 10 750 K.
In the text

	Fig. 2 SPHERE IFS YJH band image of HIP 81208 from (a) the 2019 epoch and (b) the 2022 epoch, reduced using a KLIP-based pipeline. The image here is the median taken over the entire wavelength band and shows a bright companion, HIP 81208 B, detected at 0.325 ± 0.001″ north of the primary star (averaged over the two epochs). The pixel scale is 7.46 ± 0.02 mas pixel−1 in both epochs.
In the text

Fig. 3 SPHERE IRDIS K band image of HIP 81208 from the (a) 2019 and (b) 2022 epoch observations, reduced using SpeCal in TLOCI mode. Each image here is the median taken over K1 and K2 band images of the respective epoch. Images show the companion HIP 81208 B, as well as a second companion HIP 81208 C detected at 0.325 ± 0.001″ north and 1.492 ± 0.001″ southeast of the primary star, respectively (averaged over the two epochs). The pixel scale was 12.258 ± 0.004 mas pixel−1 for K1 band and 12.253 ± 0.003 mas pixel−1 for K2 band.

In the text

	Fig. 4 K1 – K2, K1 colour–magnitude diagram, with the candidates HIP 81208B and HIP 81208C shown along with other known substellar companions, young and/or dusty objects, and field stars.
In the text

Fig. 5 Relative astrometric shifts between the two epochs along α* (x-axis) and δ (y-axis) for all the sources detected in the IRDIS FoV. The black dotted curve shows the expected motion of a pure background source with null proper motion when compared to a source co-moving with the target (the grey shaded region encompasses the uncertainty associated with this motion). The mean astrometric shifts of the probable background objects (blue crosses) detected in the data are shown as a thick blue cross, and those of HIP 81208 B and C are shown as a solid green triangle and a solid yellow square, respectively.

In the text

Fig. 6 Proper motion in αcos(δ) and δ for the sample of Gaia EDR3 sources within 5 arcmin of the primary, and for the six background sources in the IRDIS data. Also shown are the median proper motion of the six probable background objects as well as the proper motion derived for HIP 81208 B and C using the same symbols as in Fig. 5. The interesting region of proper motion defined by Eq. (2) is over-plotted in the figure as a yellow-shaded region enclosed within black dashed contours.

In the text

	Fig. 7 Contrast spectrum for HIP 81208B obtained from reduced IFS data in the YJH band and IRDIS data in K1 and K2 bands, along with the associated uncertainties. The red solid squares represent the contrast obtained from the 2022 observations and the green solid circles represent the contrast obtained from 2019 observations.
In the text

	Fig. 8 ${\chi }_{\text{red}}^2$ computed between the observed spectrum and the AMES-Cond and BT-Settl models for different effective temperatures for both 2019 and 2022 epochs. The log(ɡ) adopted for both models is 4.5. For both epochs, the models reach a ${\chi }_{\text{red}}^2\sim 1$ at Teff = 2900 K. The closest spectral type is M5.
In the text

Fig. 9 Observed spectrum for HIP 81208B in 2019 (green solid circles) and 2022 (red solid squares) epochs, along with the associated uncertainties in flux. Also plotted are the AMES-Cond 2000 log(g) = 4.5 theoretical model for Teff = 2800, 2900, and 3000 K and the BT-Settl log(ɡ) = 4.5 theoretical model for Teff = 2900 K. The best-fit models to the observed spectrum in both epochs correspond to Teff = 2900 K.

In the text

	Fig. 10 Corner plots for the accepted orbits from OFTI for companions B (a) and C (b). The median values along with upper (84%) and lower (16%) quantile intervals (~±1σ) are shown as dashed vertical lines on the posterior distributions of the respective orbital parameters. The 2D contour plots show the 1σ, 2σ, and 3σ confidence limits to these values.
In the text

	Fig. 11 Distribution of the relative inclinations between B and C’s orbits from the sample of solutions obtained by the OFTI algorithm. The light red area corresponds to systems in the ZLK configuration, or Kozai regime.
In the text

Fig. 12 Illustration of 20 random orbits drawn from the posterior distributions from the OFTI run for HIP 81208B. The left panel shows the motion of B in RA and Dec over an entire orbital period for these 20 cases. The adjacent color bar shows the epochs along the orbital positions. The yellow solid circles represent the predicted position of the companion in these orbits in the two observation epochs. The blue star represents the central primary star HIP 81208. The right panel shows the projected separation ρ and PA (deg) of B over time, as predicted from these 20 orbits. The two purple dots on both figures represent the measured values of the respective parameters at the two observation epochs.

In the text

Fig. 13 Illustration of 20 random orbits drawn from the posterior distributions from the OFTI run for HIP 81208C. The left panel shows the motion of C in RA and Dec over an entire orbital period for these 20 cases. The change in the position of the companion over the two epochs is very small and the separate data points are not visible as in the previous figure and are therefore represented by a single yellow solid circle. The right panel shows the projected separation ρ and PA (deg) of C over time as predicted from these 20 orbits. The symbols and colour bar hold the same meaning as in the previous figure.

In the text

Fig. 14 Maximum eccentricity reached by either orbit within 107 yr of evolution depending on the initial orbital elements. Black solid circles depict unstable cases (emax > 1, 44% of the cases). These results were obtained by running N-body simulations starting from 1000 solutions of the orbital fit. Configurations where orbits are crossing or almost crossing are always unstable. Orbits with relative inclination between 40 and 130 degrees are locked in a Kozai resonance, and may experience high eccentricity variations without compromising the stability of the system.

In the text