© ESO 2021

1. Introduction

Solar-type multiple systems are at least as common as individual stars: the fraction of triple-star systems was found to be 8 ± 1%, and it drops to 3 ± 1% for higher-multiplicity systems (Raghavan et al. 2010). Similarly, observations of F and G stars within 67 pc of the Sun (Tokovinin 2014a,b) show that ≈10% of all stellar systems are triple and ≈4% are quadruple. The high-order multiplicity fraction increases with stellar mass (Duchêne & Kraus 2013). Multiple star systems with n > 2 are nearly always hierarchical, meaning that they can decompose into binary or single sub-systems based on their relative separations (e.g., two close binaries that orbit each other with a wide separation). A hierarchical system can have many distinct configurations. For instance, quadruple systems can have two possible configurations. A triple system orbited by a distant fourth companion corresponds to the 3+1 configuration. The 2+2 configuration consists in two close binaries orbiting around each other. The 2+2 configuration seems to be ∼4 times more frequent than the 3+1 configuration for solar-type stars (Tokovinin 2014b). The orbital parameters in hierarchical systems could provide additional information about their formation history. It is expected that different formation processes, such as core fragmentation, disk instability, dynamical interactions, or a combination of different formation channels, leave imprints on the mass ratio, periods, eccentricities, and mutual orbit inclination of hierarchical systems (Kozai 1962; Lidov 1962; Whitworth 2001; Sterzik & Tokovinin 2002; Lee et al. 2019; Tokovinin & Moe 2020). In the last decades, observational and theoretical efforts have led to a better understanding of the formation and dynamical stability of such multiple systems (Kiselev & Kiyaeva 1980; Tokovinin et al. 2006; Eggleton 2009; Tokovinin 2018a; Hamers et al. 2021).

Wide binaries show a strong preference to be in hierarchical systems in low density young associations (Elliott et al. 2016; Elliott & Bayo 2016) and star-forming regions (Joncour et al. 2017). The fact that this relation is not seen in the same proportion in denser environments or systems in the field suggests that this could be the result of dynamical processing or the unfolding of hierarchical systems (Sterzik & Tokovinin 2002; Reipurth & Mikkola 2012). In that regard, characterising young (1–100 Myr), hierarchical systems helps to observe their early evolution. The formation channels of hierarchical systems cannot be easily determined by only characterising field stars, where billions of years of dynamical evolution may have erased their formation history. Consequently, the study of young (1–100 Myr) hierarchical systems is an important step to better understand their formation pathway. Large-scale surveys provide crucial information that helps to discover such multiple systems, but they are not well suited to finely constrain their orbital architecture. We need high-precision astrometry and radial velocity (RV) follow-up observations of the identified hierarchical systems to accurately constrain parameters of their inner and outer orbits.

The HD 98800 is a well-known hierarchical quadruple star system, and a member of the 10-Myr old TW Hydrae association (Torres et al. 2008). Located at 44.9 ± 4.6 pc from Earth according to the latest reduction of HIPPARCOS data¹ (van Leeuwen 2007), corresponding to a parallax of 22.27 ± 2.32 mas, it consists of two pairs of spectroscopic binaries (hereafter, AaAb and BaBb, see Fig. 1). Both binaries orbit each other with a semi-major axis of ≈45 au (Tokovinin et al. 2014). The AaAb system is a single-lined spectroscopic binary (SB1) with a period of 262 days (Torres et al. 1995). The mass of the Aa was estimated from pre-main sequence evolutionary models as 1.1 ± 0.1 M_⊙ (Prato et al. 2001). The BaBb subsystem is a double-lined spectroscopic binary (SB2) with a period of 315 days, the astrometric orbital solution of this binary was first presented in Boden et al. (2005) using five Keck Interferometer (KI) epochs combined with Hubble Space Telescope astrometry, and available RV observations. From this orbital solution, Boden et al. (2005) estimated a parallax of 23.7 ± 2.6 mas, and dynamical masses for Ba and Bb of 0.699 ± 0.064 and 0.582 ± 0.051 M_⊙, respectively. The BaBb pair also harbours a bright circumbinary protoplanetary disk (Skinner et al. 1992; Zuckerman & Becklin 1993), and ALMA observations revealed that the disk and the binary orbital planes are perpendicular to each other (Kennedy et al. 2019). Numerical simulations suggest that this ‘exotic’ (yet stable) configuration can be reached in some multiple systems, the so-called polar configuration (Verrier & Evans 2008; Farago & Laskar 2010; Aly et al. 2018). Dynamical evolution studies show that an inclined circumbinary disk around a highly eccentric (e ≳ 0.7) inner binary can evolve towards this configuration (Aly et al. 2015; Zanazzi & Lai 2017; Cuello & Giuppone 2019).

Fig. 1. Schematic view of HD 98800 orbital configuration. The BaBb subsystem hosts a circumbinary disk in polar configuration and is orbiting around the AaAb binary in a highly inclined orbit with a semi-major axis of ≈45 au.

Recently, the orbital characterisation of hierarchical systems hosting disks has provided new insights on the mechanism involved in the formation of multiple systems and their interaction with the disk (Kraus et al. 2020; Czekala et al. 2021). To better understand the source of disk misalignment and the formation process behind hierarchical systems, better information on well-characterised multiple systems’ architectures will be necessary. In that regard, the full characterisation of the HD 98800 quadruple system presents an opportunity to expand the sample of hierarchical systems hosting a protoplanetary disk.

In this work, we present new long-baseline infrared interferometric observations of both AaAb and BaBb subsystems, as well as new RV measurements from original observations and archival reduced spectra. The new interferometric observations resolve the relative position of Ab with respect to Aa for the first time, providing one of the missing keys for the full characterisation of this quadruple system. Additionally, we also present two new astrometric positions for BaBb, allowing us to refine the orbital solution reported in Boden et al. (2005). With the new orbital solutions of AaAb and BaBb, we re-estimated the orbital parameters of the AB outer orbit, evaluate the dynamical stability of this system, and discuss possible formation scenarios for this 2+2 quadruple.

2. Observations, astrometry, and RV

2.1. PIONIER observations and data reduction

We used the Very Large Telescope Interferometer (VLTI, Haguenauer et al. 2008; Haubois et al. 2020) with the four-telescope combiner PIONIER in the H band (1.5 − 1.8 μm, Le Bouquin et al. 2011) to observe the HD 98800 quadruple system. Our observations were carried out using the 1.8 m Auxiliary Telescopes with small and medium configurations, providing six projected baselines per configuration ranging from ∼20 to 100 m. This configuration provides an angular resolution of ∼4 mas. The estimated interferometric field of view for PIONIER is ∼160 mas (Hummel et al. 2016), but given the loss of coherence caused by spectral smearing of the companion, with our given configuration, we have a field of view ≲60 mas (Le Bouquin & Absil 2012; Gallenne et al. 2015).

The first observations of both sub-systems were taken in April and May 2019 as a part of the science verification (SV) campaign² of the New Adaptive Optics Module for Interferometry (NAOMI, Woillez et al. 2019). These observations showcase the improvement provided by NAOMI on the sharpness of the point spread function (PSF) (despite ∼1″ seeing conditions), which led to a better injection of the light in the fibre, and allowed us to mitigate light-contamination effects between A and B subsystems (A–B separation ≲0.4″). After the SV run, we obtained six more PIONIER epochs for the AaAb binary between February 2020 and March 2021 (see Table 1).

To monitor the instrumental and atmospheric transfer functions, the standard observing procedure is to interleave science and reference stars (CAL-SCI-CAL-SCI-CAL sequence). The calibrators, listed in Table A.1, were selected using the SearchCal software (Bonneau et al. 2006, 2011; Chelli et al. 2016) provided by the Jean-Marie Mariotti Center (JMMC³). The data were reduced with the pndrs package described in Le Bouquin et al. (2011). The main procedure is to compute squared visibilities (V²) and triple products for each baseline and spectral channel, and to correct for photon and readout noise. The calibrated data are available in the Optical Interferometry DataBase⁴. In Fig. 2, an example of the squared visibilities and closure phases (CP) for one of our observations of AaAb is presented.

Fig. 2. Squared visibility and closure phase measurements from one observation of AaAb taken in March 2021. The data are in blue, while the red dots represent the best binary model fitted with CANDID for this epoch. Bottom panels: residuals in the number of sigmas.

2.2. Determining the AaAb and BaBb astrometry

For each PIONIER observation, we determined the astrometric positions by fitting the V² and CP with a binary model using the interferometric tool CANDID⁵ (Gallenne et al. 2015). For each epoch, the tool delivered the binary parameters, namely the flux ratio (f₂/f₁) and the relative astrometric position (Δα, Δδ). CANDID can also fit the angular diameter of both components, however, in our case, we kept them fixed at 0.3 mas during the fitting process as the VLTI baselines did not allow us to resolve such small diameters. Briefly, the tool provides a 2D-grid of a multi-parameter fit using a least-squares algorithm (see Fig. 3). Given the small separation between AaAb and BaBb (≲0.4″), we also fitted an additional parameter to take the background cross-contamination into account, the non-coherent light, parametrised in CANDID as a resolved flux (f_res). The final astrometric positions for all epochs of each subsystem are listed in Table 1. CANDID estimates the uncertainties using a bootstrapping approach (with replacement) using 10 000 bootstrap samples. For the flux ratio and resolved flux, we used the bootstrap sample distributions and took the median value as the best-fit result and the maximum value between the 16th and 84th percentiles as uncertainty. For the astrometry, the 1σ error region of each position is defined with an error ellipse parametrised with the semi-major axis σ_maj, the semi-minor axis σ_min, and the position angle σ_PA measured from north to east. We also quadratically added the systematic uncertainty of 0.35% from the precision of the PIONIER wavelength calibration to σ_maj and σ_min (Kervella et al. 2017; Gallenne et al. 2018).

Fig. 3. Detection level map from CANDID for the observation of AaAb taken in April 2019. The colourbar shows the significance of the companion detection in the number of sigmas. The red cross points to the best-fit position.

2.3. AB astrometry

We gathered astrometric measurements from the Washington Double Star catalogue (WDS, Mason et al. 2001). The AB pair has been observed since 1909, and observations before 1991 have no reported uncertainties. For those observations, the expected astrometric uncertainty was found to be between 0.02 − 0.1″, depending on the target brightness and observing conditions (Douglass et al. 1992; Torres et al. 1999). Since 2009, the pair has been regularly observed with the speckle camera (HRCam, Tokovinin 2018b) mounted on the 4.1 m Southern Astrophysical Research Telescope (SOAR); the last observation presented in this work was obtained in April 2021.

2.4. CTIO spectroscopy

Five observations were taken with the 1.5 m telescope located at the Cerro Tololo Inter-American Observatory (CTIO) in Chile, and operated by the Small and Moderate Aperture Research Telescopes System (SMARTS) Consortium⁶, from April-July 2021. Observations were made with the CHIRON optical echelle spectrograph (Tokovinin et al. 2013). The RVs were determined from the cross-correlation function (CCF) of echelle orders with the binary mask based on the solar spectrum, as detailed in Tokovinin (2016). From these observations, we obtained five RV measurements for Aa (brighter component). The Ba and Bb components were totally blended with Aa at two epochs and they could not be separated by a multi-component fit. However, the blending certainly biases the RVs of Aa, increasing the uncertainty of these measurements. In three observations, we were able to obtain reliable RV measurements for Ba and Bb; in one of them, however, the components were still partially blended, so a larger uncertainty was assigned to it.

2.5. Reduced spectra from public archives

We found nine science-ready datasets in the ESO Phase 3 public archive⁷ taken with the Fibre-fed Extended Range Échelle Spectrograph (FEROS/2.2 m, Kaufer et al. 1999). The 1D Phase 3 spectra are given in the barycentric reference frame. One observation was taken in 2015 and the remaining eight were acquired between July and August 2007. The RVs were determined by cross-correlation with the same solar-type binary mask as used in CHIRON. The lines are blended and dominated by the lines of Aa. Consequently, with these FEROS spectra, we obtained only RV measurements for the Aa component, potentially biased by blending with Ba and Bb. Additionally, we found two reduced spectra in the ELODIE public archive⁸ at the Observatoire de Haute-Provence (OHP, Moultaka et al. 2004). The observations were taken in 1998, on January 28 and 29. The spectra are not given in the barycentre reference frame; a correction was therefore applied after retrieving the data. The RVs were determined from the CCF of the spectra with a CORAVEL-type G2 numerical mask using a standalone CCF tool⁹ (for further details, see Zúñiga-Fernández et al. 2021). The spectra from ELODIE show partially blended lines and were fitted by three Gaussian components. The RV measurements for BaBb from ELODIE have large uncertainties, but still allowed us to compute the systemic velocity of BaBb at this epoch.

2.6. Literature data

From the literature, we collected a diverse dataset for this system. The RV measurements of the primary star of AaAb (single-line spectroscopic system, SB1) and both components for BaBb (double-line spectroscopic system, SB2) were taken from Torres et al. (1995), hereafter TO95. For the BaBb binary, we also retrieved interferometric V² measurements, obtained with the KI, and published in Boden et al. (2005). Additionally, assuming that the RV of B is the same as the systemic velocity of the disk, we include the disk RV derived from the CO modelling presented in Kennedy et al. (2019), hereafter KE19.

3. Orbital fitting

We modelled our dataset with the exoplanet software package (Foreman-Mackey et al. 2020), which extends the PyMC3 framework (Salvatier et al. 2016) to support many of the custom functions and distributions required when fitting orbital parameters. Some of the parameters describing the primary or the secondary star orbits around the centre of mass are identical for both components, for example, the period (P), eccentricity (e), inclination (i), and longitude of the ascending node (Ω). But others depend on the component used as a reference, for example, the semi-amplitude of the RV (K_primary and K_secondary) and the argument of the periastron (ω_primary = ω_secondary + π). In an astrometric-only orbital fitting, it is common practice to report ω = ω_secondary, whereas with an RV orbit it is generally common practice to report ω = ω_primary. Then, in a joint astrometric-RV orbit, there could be ambiguity regarding the convention used for ω. We adopted the orbital convention from exoplanet¹⁰, where the argument of periastron ω is reported with respect to the primary star, and the longitude of the ascending node Ω is the node where the secondary is moving away from the observer (see Fig. 4).

Fig. 4. Diagram of the orbit of the secondary star around the centre of mass (yellow plane) and the reference plane (grey). This diagram follows the orbital convention of exoplanet.

Given that BaBb is an SB2, the orbital fitting procedure is slightly different compared to the AaAb subsystem (SB1). The different steps for each orbital fitting are explained below. The Markov chain Monte Carlo (MCMC) samples and PyMC3 models corresponding to both subsystems are available online¹¹. The prior distributions and corner plots from the orbital parameters’ posterior samples are displayed in Appendix B.

3.1. BaBb orbit

Given that BaBb is an SB2, we can fit the astrometric points from PIONIER together with the RV amplitude of each component, K_Ba for the primary and K_Bb for the secondary, as well as the systemic velocity γ_B. Additionally, we extended exoplanet to include the V² model for individually unresolved components in a binary system. Briefly, the fringe contrast V² of a binary system depends on the properties of the individual components and the binary separation (Berger & Segransan 2007),

$$\begin{array}{c}\hfill V_{\mathrm{binary}}^2=\frac{1+{\left(\frac{f_2}{f_1}\right)}^2+2\left(\frac{f_2}{f_1}\right)cos\left(\frac{2\pi C(u\mathrm{\Delta }\alpha +v\mathrm{\Delta }\delta )}{\lambda }\right)}{{(1+\frac{f_2}{f_1})}^2},\end{array}$$(1)

where Δα and Δδ are the relative separation in right ascension and declination, respectively (from the exoplanet model), u and v are the projected baselines (in meters), f₂/f₁ is the flux ratio, and λ is the wavelength. The parameter C is a conversion factor so that the astrometry is in arcsec and the wavelength in μm. The V² measurements from KI were then included in the fitting process, where the flux ratio f₂/f₁ was fitted as a free parameter along with the other orbital parameters (see Table 2).

All the orbital parameters were estimated from the posterior distributions, taking the median values as the best-fit results and the maximum values between the 16th and 84th percentile as uncertainties. From these distributions, we could then calculate the distribution of the masses for both components as well as the distance to the system using the following equations (Torres et al. 2010; Gallenne et al. 2019):

$$\begin{array}{cc}\hfill M_{\mathrm{Ba}}& ={\displaystyle \frac{1.036149\times {10}^{-7}{(K_{\mathrm{Ba}}+K_{\mathrm{Bb}})}^2{K}_{\mathrm{Bb}}P{(1-e^2)}^{3/2}}{{sin}^{3}i}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill M_{\mathrm{Bb}}& ={\displaystyle \frac{1.036149\times {10}^{-7}{(K_{\mathrm{Ba}}+K_{\mathrm{Bb}})}^2{K}_{\mathrm{Ba}}P{(1-e^2)}^{3/2}}{{sin}^{3}i}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill a_{\mathrm{au}}& ={\displaystyle \frac{9.191940\times {10}^{-5}(K_{\mathrm{Ba}}+K_{\mathrm{Bb}})P\sqrt{1-e^2}}{sini}},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \pi & ={\displaystyle \frac{a}{a_{\mathrm{AU}}}},\hfill \end{array}$$(5)

where M_Ba and M_Bb correspond to the masses of the primary and the secondary stars, respectively, expressed in solar mass, P is the period in days, K_Ba and K_Bb are the RV semi-amplitudes of the primary and secondary star in km s⁻¹, respectively, and a is the angular semi-major axis in arcseconds. The parameter a_au is the semi-major axis expressed in astronomical units. Table 2 lists a full description of the inferred orbital parameters. Figure 5 shows the best-fit RV curve. Figure 6 shows the best-fit visual orbit; the black dots are the phase coverage of the KI observations, that is the astrometric positions from the best-fit orbit corresponding to the observation date of each V² dataset (see Fig. B.1). Some parameters seem incompatible with the previous result taking the uncertainties into account; this may be due to the fact that some of the uncertainties could have been underestimated.

Fig. 5. Phase-folded RVs orbit for BaBb. The systemic velocity γ for each set of observations was subtracted. The solid line corresponds to the best-fit model. Upper panel: RVs of Ba, and lower panel: Bb.

Fig. 6. Best orbital solution for BaBb. The solid line corresponds to the best-fit model and the shaded area to the 1σ region. The primary star Ba is located at the origin. The relative positions of Bb are plotted as filled dots. The error ellipses from PIONIER astrometry are smaller than the markers.

3.2. AaAb orbit

This subsystem is an SB1, therefore it is not possible to break the degeneracy between the parallax and the semi-major axis and determine individual stellar masses. The orbit is based on the astrometric points from PIONIER and on the RVs of the primary star Aa. Consequently, we only fitted the RV semi-amplitude of the Aa component of the system K_Aa, and the systemic velocity γ_A. To estimate the masses of the individual components of an SB1, we must assume a distance. We tested two parallax values, the one obtained from the orbital fitting of BaBb and the HIPPARCOS one (van Leeuwen 2007), as there is currently no reliable Gaia parallax for the system. In our MCMC model, we included these parallax values as a prior using a normal distribution (22.27 ± 2.31 mas and 22.0 ± 0.6 for the HIPPARCOS one and the one based from the orbital solution of BaBb, respectively). The parallax is then a free parameter in our orbital fitting using the abovementioned priors. Using Kepler’s third law and combining Eqs. (3) and (4), we calculated

$$\begin{array}{cc}\hfill M_{\mathrm{tot}}& ={\displaystyle \frac{a_{\mathrm{AU}}^3}{P_{\mathrm{years}}^2}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill M_{\mathrm{Ab}}& ={\displaystyle \frac{1.036149\times {10}^{-7}{K}_{\mathrm{Aa}}\sqrt{1-e^2}{a}_{\mathrm{AU}}^2}{{(9.191940\times {10}^{-5})}^{2}Psini}},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill M_{\mathrm{Aa}}& =M_{\mathrm{tot}}-M_{\mathrm{Ab}},\hfill \end{array}$$(8)

where M_tot, M_Aa, and M_Ab correspond to the total mass, and the primary and secondary star masses, respectively, expressed in solar mass, P_years the period in years, P the period in days, K_Aa the RV semi-amplitude of the primary star in km s⁻¹, and a_AU the semi-major axis expressed in astronomical units.

The posterior distributions for the masses and parallaxes, assuming different initial priors for the parallaxes are shown in Fig. 7 in blue and red, respectively. The orbital parameters converge and have the same results for both cases, only the physical parameters that are dependent on the distance are affected by the choice of the prior distribution for the parallax (i.e., M_Aa, M_Ab, and a_AU).

Fig. 7. Posterior distribution of masses and parallax of AaAb subsystem assuming the HIPPARCOS value (red) and the solution from BaBb fitting (blue) as a prior distribution of the parallax in our MCMC model. The red and blue lines highlight the median of each distribution.

All the orbital parameters were estimated from the posterior distributions taking the median values as the best-fit results and the maximum values between the 16th and 84th percentile as uncertainties (see Table 3). Figure 8 shows the best-fit binary orbit (identical in the plane of the sky for both parallax scenarios). In the rest of the paper, we assume the masses of AaAb as derived with the parallax obtained from the BaBb best orbital fit.

Fig. 8. Best orbital solution for AaAb. In both panels, the solid line corresponds to the best-fit model. Bottom panel: the primary star Aa is located at the origin. The relative positions of Ab are plotted as filled dots; the error ellipses from PIONIER astrometry are smaller than the marker. Upper panel: the coloured markers correspond to the primary star RV measurements. The systematic velocity γ for each set of observations was subtracted.

3.3. Outer orbit A–B

Using our orbital solutions of the inner binaries of the system, we recalculated the orbital parameters of AB. We assume that the systemic velocities of AaAb and BaBb from Torres et al. (1995), FEROS, CHIRON, and ELODIE observations in our orbital fitting results, and the one from CO modelling by Kennedy et al. (2019) (KE19), correspond the centre-of-mass RVs of A and B in the outer orbit (see Table A.4). We jointly fitted the astrometric position with the RV measurements of AB. In our MCMC model, we included the parallax and the masses obtained from the inner orbits’ results as priors. For consistency, we used the AaAb masses derived from the parallax obtained from the orbital fitting of BaBb. The normal distribution priors for the masses and parallax are M_A : 1.22 ± 0.5 M_⊙, M_B : 1.38 ± 0.5 M_⊙, and π : 22.0 ± 0.6 mas, respectively. The γ_AB was included as a free parameter, with a uniform prior between 0 and 20 km s⁻¹. Given that the visual micrometric measurements made before 1991 have unknown uncertainties, we defined the large (σ ∼ 0.1″) and the small (σ ∼ 0.02″) uncertainty cases for these measurements (solutions I and II in Table 4), according to the typical range of errors reported in the astrometry measurements by USNO (Douglass et al. 1992; Torres et al. 1999).

All the orbital parameters were estimated from the posterior distributions, taking the median values as the best-fit results and the maximum values between the 16th and 84th percentile as uncertainties (see Table 4). The best-fit orbits for solutions I and II are shown in Figs. 9 and 10, which are in good agreement with the astrometric measurements. We show the phase-folded RV best-fit orbit in Fig. 11; the narrow 1-sigma region in this figure mainly comes from the constraints imposed by the AB masses and parallax prior distributions. There are likely small instrumental zero-point offsets among the data sets that were used to determine the systemic velocity variation, which are difficult to determine and could bias the outer orbit solution. As a reminder, all these results rely on the masses and parallax estimated in the orbital fitting of the inner subsystems. The parallax and masses derived from BaBb RV semi-amplitudes (K_Ba and K_Bb) are proportional to K² and K³, respectively. Thus, a small systematic error in K_Ba or K_Bb can bias the masses and parallax results substantially. The RV amplitudes may be biased, especially for the weakest lines of Bb. Therefore, the masses and the parallax of the BaBb pair that mainly rely on the RVs by Torres et al. (1995) should be considered with caution. A small change in the method of splitting the blended spectra can lead to different masses. The posterior samples of the orbital parameters and all prior distributions used in the MCMC model are available in Appendix B.

Fig. 9. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. The solid black line corresponds to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue one assumes large uncertainties (σ ∼ 0.1″). The red dots correspond to the astrometric measurements. For better visualisation, the error bars are shown only in Fig. 10.

Fig. 10. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. In all panels the solid line corresponds to the best-fit model and the shaded area to the 1σ region. The solid black lines correspond to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue ones assume large uncertainties (σ ∼ 0.1″). The red dots correspond to the astrometric measurements. The error bars shown in the astrometry before 1991 correspond to the large uncertainty case.

Fig. 11. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. In both panels the solid line corresponds to the best fit model and the shaded area to the 1σ region. The solid black lines correspond to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue ones assume large uncertainties (σ ∼ 0.1″). The dots markers correspond to the RV measurement of systemic velocities from our orbital solutions and the one obtained from CO modelling.

Accurate astrometry of AB reveals a wavy motion caused by the subsystems (wobble); its amplitude gives an independent constraint on the inner mass ratios. Neglecting the smaller wobble of BaBb, we modelled the astrometry of AB by a combination of two Keplerian orbits, with the orbital parameters of AaAb fixed to the values determined above. A simple least-squares fit yielded the AB orbital parameters similar to solution I, for example P = 233 ± 41 yr. The ratio of the wobble amplitude to the semi-major axis of AaAb was found to be f = 0.18 ± 0.04. Neglecting the influence of the faint light of Ab on the photo-centre of A, this factor gives the inner mass ratio q = f/(1 − f) = 0.22, compatible within errors with the mass ratio of 0.31 estimated above from the orbit of AaAb.

4. Short- and long-term future of the quadruple system

Several studies investigated the stability of the system over time (e.g., Verrier & Evans 2008; Kennedy et al. 2019), but those studies mostly focused on the stability of the disk around BaBb and less about the evolution of the quadruple system itself. In this section we intend to study both the short- and long-term dynamical evolution of the two pairs of binary systems. Using the new (or revised) orbits obtained for AaAb, BaBb, and AB, we first quantify the dynamical stability over time of the four stars, and second make preliminary predictions for the transit of BaBb and its disk in front of AaAb. To make such predictions, we use the N-body code REBOUND¹² (Rein & Liu 2012). For the simulations, we used the orbital solutions for AaAb and BaBb listed in Tables 3 and 2, and we tested both solutions I and II for the orbit of AB (Table 4). The best-fit parameters are directly taken from the posterior distributions, as their median values.

4.1. Dynamical stability

For our dynamical stability analysis, we use the ‘mean exponential growth of nearby orbits’ (MEGNO) criterion implemented in the REBOUND package. As discussed in Hinse et al. (2010), the MEGNO factor, first introduced in Cincotta & Simó (2000), provides an estimate of how ordered or chaotic a system is. The MEGNO factor is the integral of variational vectors for a given integration time and a given set of parameters. It is therefore necessary to sample different timescales as the MEGNO is expected to vary over time (and converge to a value of 2 for a stable system), tracing the different orbital timescales. In our case, we want to study the stability of the orbits by changing the masses of Aa and Ab. To do so, we computed the MEGNO value for a matrix of masses, the rows of the matrix consist of 12 linearly spaced masses for M_Aa in the range [0.5, 1.5] M_⊙ and 10 linearly spaced masses for M_Ab in the range [0.1, 1.0] M_⊙.

To setup the simulation, we sequentially added Aa and Ab, and then included a third particle representing BaBb as a single star. We then integrated the motion of all stars forward in time, using the IAS15 integrator (Rein & Spiegel 2015). Since it is necessary to capture the different timescales for the evolution of the system, we used several integration times, in years: 1000, 2000, 5000, 10 000, 20 000, 50 000, 100 000, 250 000 500 000, and 1 000 000. For all the simulations, none of the stars escaped the system, suggesting that it is stable over thousands of AB orbits. Since the final matrices all are homogeneous, in Fig. 12 we show the mean MEGNO value (over the 12 × 10 matrix) as a function of the integration time, for both solutions I and II, and we computed the standard deviation as the uncertainties. The stability criterion displays an exponential decay, converging towards a value of 2 (Hinse et al. 2010), therefore indicating that the system should be stable over a long period of time, regardless of the uncertainties on the AB orbital parameters. In the exercise above, we treated BaBb as a single star, and it might be worth re-visiting the stability of the system when considering all four stars, but given the large uncertainties on the AB orbit, this is out of the scope of this study.

Fig. 12. Mean MEGNO value for the 12 × 10 matrix for different integration times. Error bars correspond to 1σ.

4.2. The transit of the disk in front of AaAb

The orbital parameters of AB strongly suggest that the BaBb pair and its disk will pass in front of the AaAb system (Kennedy et al. 2019), starting sometime in 2026 (depending on the solution used for AB). This presents a unique opportunity to observe and characterise the properties of the dust and gas disk via photometric (and spectroscopic) monitoring of the whole system. In virtue of this possible occultation, we investigate how the photometric light curve might look, including the four stars in the simulation to account for possible interactions between the two binary systems (in App. C we provide a more detailed explanation on how the simulation is initialised).

To make the predictions for the transit, the starting time of the simulation was set to 2015.17. The choice of the starting date does not matter as we are using the orbital solutions determined in this paper. We integrated the simulation over 18 years and saved 10 000 intermediate steps (one every ∼0.7 days), saving the positions of the four stars in the reference system centred at the centre of mass of BaBb. The integration was done using the IAS15 integrator, but we also compared our results with the WHfast integrator (Rein & Tamayo 2015) with a timestep of 0.011 days, and found no significant differences between the two simulations. Additionally, saving more intermediate steps does not lead to a significant improvement of the resolution of the simulated light curve.

With the (x, y) positions of Aa and Ab, on the plane of the sky, we then estimated if they overlap with the position of the disk, which is centred at the centre of mass of the B system (Fig. C.1). We used the parameters reported in Kennedy et al. (2019), namely, the inner and outer radii (2.5 and 4.6 au, respectively), eccentricity (0.03), position angle (15.6°), inclination (26°), and argument of periapsis (−73°). To estimate the extinction caused by the circumbinary disk, we first needed the flux ratio between Aa and Ab, and an analytical form for the integrated vertical optical depth of the disk. For the flux ratio, we used the results from the modelling of the PIONIER observations, and the normalised fluxes are F_Aa = 0.87 and F_Ab = 0.13 in the H band. For the vertical optical depth, it is parametrised as τ(r) = 0.5 × r₀/r, where r₀ is the inner radius of the disk (τ(r) = 0 inside and outside the disk). Before applying the extinction law, we first needed to estimate the distance r in the midplane of the disk, accounting for projection and rotation effects. We therefore defined a rotation matrix ℛ based on the inclination, argument of periapsis, and position angle of the disk, and de-projected the on-sky (x, y) positions of the disk and Aa and Ab stars. The normalised flux at each time-step is then F_Aae^−τ(r_Aa) + F_Abe^−τ(r_Ab) (the contribution of BaBb is neglected here, but since the vertical optical depth of the disk remains unknown the absolute depth of the transit cannot be constrained).

We then simulated 1000 transits and their respective light curves by modifying the AB orbital parameters (for both solutions I and II), the parallax, and all four masses randomly drawing 1000 realisations from the MCMC fitting of the AB orbit. This ensures that the correlations between the different parameters are preserved (to avoid, for instance, a small semi-major axis and large stellar masses that would lead to a much shorter orbital period). Finally, from these 1000 light curves, we estimated a probability distribution of the normalised flux as a function of time, which is shown in Fig. 13, where the transit light curve for the best-fit solution is shown in orange. The figure shows that the transit should be well constrained in time, and we predict it to start in 2026, going out and passing through the inner regions (devoid of dust) before re-entering behind the northern side of the disk. Our simulations suggest that the transit event should finish sometime between 2030 and 2031. The best-fit solution shows the complex structure of the light curve as one of the stars is sometimes not occulted by the disk. Comparing the light curves for both solutions I and II, we note that transit starts earlier for solution I, but both cases show a similar behaviour. Overall, regular photometric monitoring of the quadruple system between 2026 and 2031 at different wavelengths would put unique constraints on the vertical optical depth of the circumbinary disk around BaBb, offering the opportunity to directly measure the surface density of the dust and to possibly derive constraints on the typical size of the dust particles.

Fig. 13. Probability density plot of 1000 realizations of the light curve for the occultation of AaAb behind the disk surrounding BaBb for solutions I and II (top and bottom, respectively). The colour bar shows the probability of getting a determined flux at a given time, such that the sum along each of the columns is normalised to unity. In orange we show the light curve for the best-fit parameters (Table 4).

5. Discussion

Here we discuss the implications of our results in the context of the formation of this quadruple system and its influence on the disk evolution. A further dynamical simulation of this system is beyond the scope of this paper.

5.1. Comparisons with previous results

We refined the orbital results from Boden et al. (2005) and resolved the orbit of the AaAb subsystem for the first time using PIONIER observations (see Tables 2 and 3). Using our orbital solution of BaBb, we derived a dynamical parallax of 22.0 ± 0.6 mas corresponding to a distance of 45 ± 1 pc. Boden et al. (2005) placed the system at 42.2 ± 4.7 pc using their orbital solution, and the updated reduction of the HIPPARCOS data (van Leeuwen 2007) measured a parallax of 22 ± 2 mas, corresponding to a distance of 45 ± 5 pc likely biased by the unresolved A–B components. There are two entries at the Gaia EDR3 catalogue (Gaia Collaboration 2021) at ∼0.1″ and ∼0.3″ from the positions of AaAb and BaBb, respectively, corresponding to the angular distance after correction using Gaia EDR3 proper motion from J2000 to J2016. Additionally, both subsystems were identified in the cross-matched catalogue between Gaia EDR3 and the Tycho-2 merged with the TDSC (I/350/tyc2tdsc, Marrese et al. 2021). The parallax values from Gaia EDR3 are 20.1 ± 0.3 mas and 23.7 ± 0.4 mas for BaBb and AaAb, respectively. However, both measurements have a large re-normalised unit weight error (RUWE) value Lindegren et al. (2018) and then are considered unreliable. The RUWE value is expected to be around 1.0 for a good fit to the astrometric observations, while in this case it is ∼9 and ∼6 for AaAb and BaBb, respectively, meaning that in both cases the unresolved companions produce motion in the photo-centre, so the 5-parameter Gaia astrometric model performs poorly. The distance inferred with our new results remains consistent with Boden et al. (2005) within 2.3σ and is compatible with the HIPPARCOS value. Using the new distance of BaBb and the orbital solution of AaAb, we derived, for the first time, the dynamical masses of the AaAb binary as M_Aa = 0.93 ± 0.09 M_⊙ and M_Ab = 0.29 ± 0.02 M_⊙. Using the Baraffe et al. (2015) 10 Myr isochrones and the dynamical masses of AaAb, we estimated an H-band flux ratio of 15.85% and 15.06% for solar and sub-solar ([M/H] = −0.5) metallicity, respectively (see Appendix D). These flux ratio values are compatible with the flux ratio derived with our PIONIER observations (∼14%, see Table 1).

Prato et al. (2001) compared the stellar properties derived from near- and mid-infrared diffraction-limited imaging with pre-main sequence evolutionary tracks, yielding masses of M_Aa = 1.1 ± 0.1 M_⊙, M_Ba = 0.93 ± 0.08 M_⊙, and M_Bb = 0.64 ± 0.1 M_⊙ and an age of ∼10 Myr. These values are compatible with the dynamical masses derived in this paper within ∼1.5σ. On the other hand, the SED models presented in Boden et al. (2005) suggested stellar properties compatible with the ones published in Prato et al. (2001). However, the predicted masses from evolutionary tracks were significantly higher than the dynamical masses from Boden et al. (2005). The authors claimed that this discrepancy came from the assumption of solar abundances in the evolutionary models, proposing sub-solar abundances ([M/H] = −0.5) with an age in the range 8 − 20 Myr. Later, Laskar et al. (2009) estimated a metallicity of [M/H] = −0.2 ± 0.1 using high-resolution echelle spectra. Additionally, they determined the visible-band flux ratio for Bb/Ba to be 0.416 ± 0.005. This value is compatible with the visible-band flux ratios estimated from Baraffe et al. (2015) isochrones at 10 Myr and our BaBb dynamical masses results of 0.458 and 0.428 for solar and sub-solar ([M/H] = −0.5) metallicity, respectively (see Appendix D). Given the uncertainty on the derived dynamical masses due to the degeneracy between K_Ba and K_Bb, we cannot use the quadruple system yet to benchmark evolutionary track models, calling for additional observations to better constrain both the orbital solutions and the abundances of the four stars. Both I and II AB orbital solutions feature comparable values for the inclination and longitude of the ascending node Ω, within ≲0.5° from the latest orbital solution (Kennedy et al. 2019, see Table 4). This result shows that despite the fact that the orbit of AB will remain uncertain for several years as more observations are collected, its orientation is already well constrained and robust.

5.2. Mutual alignment

The mutual inclinations between the inner and outer orbits in a hierarchical system can constrain the initial conditions of its formation (Fekel 1981; Sterzik & Tokovinin 2002). Hierarchical fragmentation of a rotating cloud (Bodenheimer 1978) or fragmentation of a circumbinary disk (Bonnell & Bate 1994) should result in near co-planar configurations. On the other hand, misaligned orbits could be the result of turbulent fragmentation or dynamical interactions (Lee et al. 2019). Similarly, the relation between circumbinary disk orientation and the orbital parameters of the host binary can be used to better constrain their formation scenarios (Czekala et al. 2019). The relative inclination Φ between the inner and outer orbits (or disk) is given by

$$\begin{array}{c}\hfill cos\mathrm{\Phi }=cos{i}_{1}cos{i}_2+sin{i}_{1}sin{i}_{2}cos({\mathrm{\Omega }}_1-{\mathrm{\Omega }}_2),\end{array}$$(9)

where i₁,  i₂ are the inclinations of each orbit (or disk and orbit) and Ω₁,  Ω₂ are the corresponding longitudes of the ascending nodes. The mutual inclination Φ ranges from 0° to 180°, where Φ = 0° corresponds to co-planar and co-rotating orbits. When Φ > 90° the systems are retrograde, and Φ = 90° means polar configuration. The circumbinary disk of BaBb was initially thought to be co-planar with the host binary (Tokovinin 1999; Prato et al. 2001), but recent ALMA observations revealed that the circumbinary disk is actually in polar configuration (Kennedy et al. 2019). Additionally, Giuppone & Cuello (2019) suggested that the near polar configuration between the circumbinary disk and BaBb orbit is the most stable configuration among all possible disk inclinations. Given that we reduced the uncertainty of i and Ω for the BaBb orbit from ∼3° to ∼0.5°, it is important to re-calculate the mutual inclination. Kennedy et al. (2019) found that the disk is inclined either by 26° or 154° with respect to the sky plane. The Ω_disk published in Kennedy et al. (2019) is defined as the node where rotation of the disk is moving towards the observer, that is with a difference of π with respect to our convention. For consistency, we added π to the published value resulting in Ω_disk = 196° ±1°. The new mutual inclinations between all the orbital planes of HD 98800 are reported in Table 5, including the mutual inclination of the disk with respect to the inner and outer binaries. For these mutual inclination values, we used the AB orbital parameters from solution I (see Table 4). Given that the inclination and Ω value of solutions I and II are close to each other within ∼0.1°, the subsequent analysis remains valid for both outer orbit solutions. The uncertainties were calculated using a Monte-Carlo uncertainty propagation, assuming Gaussian errors. We confirmed the near polar configuration of the disk relative to the orbital plane of BaBb in the case i_disk = 26° and found Φ_{BaBb − Disk} = 134.2° in the case i_disk = 154°. Using the posterior distributions from our fitting, and the posteriors from the disk fitting from Kennedy et al. (2019), yields a nominally significant misalignment of the disk angular momentum vector and the binary pericentre vector; 1.7 ± 0.5°. In principle, this misalignment provides a measurement of the disk mass, but given likely systematic uncertainties, for example, in estimating blended RVs, we consider this measurement to be an upper limit. Updating the calculation from Martin & Lubow (2019) using the 99.7th percentile from our posteriors, the upper limit on the disk mass is 0.02 M_⊙. The angle from polar is slightly smaller, but the binary mass is larger, so our limit is essentially the same as the upper limit computed by Martin & Lubow. Circumbinary disks are preferentially co-planar around short period (< 40 days) host binaries (Czekala et al. 2019), while for longer orbital periods, mutual inclinations are found in a wide range of configurations (Kennedy et al. 2019; GRAVITY Collaboration 2021; Czekala et al. 2021). In that regard, determining the orbital parameters of binaries and the mutual inclination of the circumbinary disk at intermediate periods (40–300 days), such as the presented HD 98800 results, can contribute to better understand the dynamical scenario leading to co-planar or polar disk configurations.

5.3. Formation history of HD 98800

The HD 98800 system is a member of the TWA Hydrae young loose association (Torres et al. 2008), therefore it is unlikely that it experienced strong external dynamical interactions with other stars. In general, hierarchical systems that formed under high dynamical interactions between nascent protostars have misaligned and eccentric orbits, and their masses are not comparable (Sterzik & Tokovinin 2002). The AB and AaAb orbits are moderately misaligned (see Table 5) and, excluding the Ab component, the masses are comparable. It is expected that the orbital and physical parameters of this quadruple system contain imprints of its formation scenario. Near co-planarity and comparable masses in wide solar-type hierarchical systems can be a sign of their formation from a common core (Tokovinin 2020a). The collapse of two nearby clouds and their inward accretion-driven migration by accretion (Tokovinin & Moe 2020) can result in compact co-planar hierarchical systems with moderate eccentricities and period ratios. However, HD 98800 is a quadruple system with a 2+2 configuration where the inner orbits are counter-rotating and the BaBb is misaligned with the outer orbit AB.

The encounter of two clumps can create shock fronts that lead to the fragmentation of each core into a binary, forming a 2+2 quadruple system (Whitworth 2001). Hypothetically, this formation scenario can produce wide quadruple systems with similar masses between all four components and comparable inner periods, called ϵ Lyr type (Tokovinin 2008), where the inner orbits are expected to be mutually misaligned. Generally, ϵ Lyr type have wide outer separations (P_outer ≳ 450 kyr), but more compact 2+2 systems are known as well (HIP 41171, P_outer ∼ 900 yr, Tokovinin 2019: FIN 332, P_outer ∼ 3000 yr, Tokovinin 2020b). Although the outer period of HD 98800 is shorter than usual for these systems (P ∼ 200 − 400 yr), the orbital configuration still matches this ϵ Lyr type except for the expectation of similar masses of its components. The mass-ratio of BaBb and AB are ∼0.8 and ∼0.9, respectively, while the mass-ratio of AaAb binary is ∼0.31.

The large BaBb eccentricity and its counter-rotating configuration with respect to the AaAb and AB orbits could be explained as the result of dynamical interactions. Tidal forces may have ripped away circumbinary material from AaAb, and in the same way, may have perturbed the BaBb circumbinary disk and the eccentricity of the host binary. In consequence, the formation process of the HD 98800 system remains unclear.

5.4. The low mass ratio of AaAb and its lack of a disk

An intriguing characteristic of HD 98800 is that it still holds a massive circumbinary gas disk around the system BaBb (Ribas et al. 2018; Kennedy et al. 2019). Nonetheless, no circumbinary disk has been found around the system AaAb. A possible explanation for the persistent existence of the detected disk has been proposed by Ribas et al. (2018). These authors speculated that the disk has survived for so long because of the tidal torques exerted by BaBb on the inner edge and by AB on the outer edge, which stopped or significantly reduced viscous accretion, leading to a scenario in which the disk is only losing mass due to photo-evaporation. On the other hand, the lack of a disk around system A, which could have evolved in a similar way as the disk around B, could be related to a faster disk dispersal due to a higher X-ray luminosity, estimated to be ∼4 times the one of system B (Kastner et al. 2004).

Recently, with a 1D+1D model of gas disk evolution, Ronco et al. (2021) explored the scenario proposed by Ribas et al. (2018) in arbitrary hierarchical triple star systems and, particularly, in HD 98800. They show that the current age and mass of gas of HD 98800 B can be reproduced if the disk was originally an intermediate to high-mass disk (∼0.05 − 0.1 M_⊙), and if it had a moderate to low viscosity (10⁻⁴ − 10⁻³). To evaluate the current non-existence of a disk around system A, these authors considered, for simplicity, that both the disk parameters and the characteristics of system A (i.e., its mass ratio and separation) were the same as those of the system B, except for its higher X-ray luminosity, as suggested by Kastner et al. (2004). Under these considerations, their simulations show that the possible disk around A may have dissipated in less than 7 − 10 Myr, the estimated age of HD 98800. We know that the assumption of equal inner mass ratios in HD 98800 does not hold. However, Ronco et al. (2021) also show that the smaller the mass ratio of the inner binary in a hierarchical triple star system, the faster the circumbinary disk dissipates, suggesting that the disk around system A in HD 98800 may have dissipated even faster. Our new findings and the characterisation of system A, presented in Sect. 3.2, effectively show a mass ratio that is much lower than that of system B, reaffirming this possibility and contributing to the explanation of the absence of the A disk.

6. Summary and conclusions

In this work, we present a new orbital solution for the HD 98800 quadruple system. Using PIONIER observations, we obtained new astrometric positions and a flux ratio of AaAb and BaBb subsystems. We refined the orbital solution presented by Boden et al. (2005) and derived, for the first time, the full orbital solution for the AaAb binary. From our orbital solution, we confirmed the polar configuration of the circumbinary disk around BaBb. Using the dynamical parallax of BaBb, we calculated the dynamical masses of the AaAb pair. The dynamical masses and parallax are strongly dependent on the RV semi-amplitude K_Ba and K_Bb, estimated mainly from the RV measurements by Torres et al. (1995). New high-resolution spectroscopic observations of HD 98800 could remove possible biases in the estimation of the RV semi-amplitude of the inner systems. Spectroscopic observation with adaptive optics correction could allow us to acquire resolved spectra of each subsystem, thus avoiding line blending of the four components. The estimated visible-band AaAb flux ratio is ≲1% (Laskar et al. 2009), making it difficult to disentangle the RV of Ab. From our PIONIER observations, we estimated an H-band flux ratio of ∼14% for the AaAb binary. This more favourable flux ratio opens the possibility to measure, for the first time, the RV of Ab using high resolution infrared spectroscopy. This would allow us to calculate the dynamical masses and parallax of Aa and Ab independently from the parallax of BaBb. Spectroscopic monitoring of HD 98800 is relevant not only for more robust dynamical masses and parallax estimates, but also to properly establish the abundances of the four stars. These measurements will provide valuable inputs to test and improve pre-main sequence evolutionary models and better constrain models of dust disk evolution.

We tested the dynamical stability of the quadruple using N-body simulations. Using the orbital parameters and the mass values of the inner binaries, the simulation probed the long-term stability of this system for both outer orbit solutions; we found that the system should be stable over thousands of orbital periods. The AB outer orbit predicts that the AaAb binary will pass behind the disk around BaBb in the coming years. Using our N-body simulation, we predicted that the transit will start in 2026 and should finish between 2030-2031. This transit presents an unprecedented opportunity to characterise the disk structure along a ∼10 au long chord, with the width of this chord set by the projected extent of the AaAb orbit.

From mass ratios, periods, eccentricities, and mutual orbit orientations, we evaluate possible formation scenarios for HD 98800. The similarity of the components’ masses suggests a common formation history. The misalignment between the orbital planes of the inner binaries and the high eccentricity of the BaBb pair suggest a possible dynamical perturbation. Assuming AaAb as a binary exactly equal to that of system B, but with a higher X-ray luminosity as suggested by Kastner et al. (2004), simulations from Ronco et al. (2021) show that the disk around A can dissipate in less than 10 Myr due to photo-evaporation. This scenario can explain the lack of a circumbinary disk around the AaAb subsystem. These authors also show that a lower mass ratio could indeed promote faster photo-evaporation of the disk. Thus, the low mass ratio derived here actually agrees with faster disk dispersal.

With the current observational evidence, we cannot properly establish the formation process of HD 98800 as there are still some uncertainties in the parallax of A as well as in the orbit of AB. Recently, other works have also used long-baseline infrared interferometry to characterise hierarchical multiple systems (Kraus et al. 2020; GRAVITY Collaboration 2021; Czekala et al. 2021). Further monitoring of other hierarchical systems, especially at young ages (1–100 Myr), in combination with large survey data, will improve our understanding of the formation and dynamical evolution of these kinds of systems.

Acknowledgments

The authors would like to thank the anonymous referee for constructive comments that helped to improve the content and clarity of this paper. S.Z.-F. acknowledges financial support from the European Southern Observatory via its studentship program and ANID via PFCHA/Doctorado Nacional/2018-21181044. J.O. acknowledges support from the Universidad de Valparaíso and from Fondecyt (grant 1180395). S.Z.-F., A.B., J.O. and M.P.R. acknowledge support by ANID, – Millennium Science Initiative Program – NCN19_171. A.B. acknowledges support from Fondecyt (grant 1190748). M.P.R. acknowledges support from Fondecyt (grant 3190336). G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. This research made use of exoplanet (Foreman-Mackey et al. 2020) and its dependencies (Astropy Collaboration 2013, 2018; Salvatier et al. 2016; Theano Development Team 2016). This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory. 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 publication makes use of VOSA (Bayo et al. 2008), developed under the Spanish Virtual Observatory project supported by the Spanish MINECO through grant AyA2017-84089. This research has made use of NASA’s Astrophysics Data System.

References

Appendix A: Observations

This section presents complementary information regarding the observations used in this work. The calibrator stars used in our PIONIER observations are listed in Table A.1. These calibrators were chosen using the SearchCal tool.

Most of the RV measurements used in this work were published by Torres et al. (1995). Here we present the new RV measurements from CHIRON observations and science ready archive spectra, see Table A.2 and Table A.3. The RVs used in the orbital fitting of the AB orbit are listed in Table A.4.

The AB astrometric measurements before 2016 are available at the Washington Double Star catalogue (WDS, Mason et al. 2001) and Tokovinin (2018c). The new astrometry measurement from speckle interferometry at SOAR are listed in Table A.5.

Appendix B: Orbital fitting complementary information

This section presents the prior distributions used for each orbital fitting. Additionally, we also show the corner plots from the posterior samples of each MCMC model. Fig. B.1 shows the V² from KI observations and the best fit binary model from the BaBb orbital fitting result.

Fig. B.1. Squared visibilities from Keck Interferometer observations published in Boden et al. (2005). The black circles represent the observed values and the red crosses represent the best-fit BaBb binary model from this work.

Fig. B.2. Posterior samples of AaAb orbital parameters. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.

Fig. B.3. Posterior samples of BaBb orbital parameters. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.

Fig. B.4. Posterior samples of AB orbital parameters for solution I. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.

Fig. B.5. Posterior samples of AB orbital parameters for solution II. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.

Appendix C: N-body simulations

In REBOUND, particles (in that case, stars) are added sequentially to the simulation. Even though a ‘primary’ keyword can be provided to indicate, for instance, that star #4 is orbiting star #3, the orbital parameters of the AB orbit are obtained with respect to the centres of mass of AaAb and BaBb, respectively. Therefore, to initialise the simulation, we determined the initial conditions of the four stars. We first added Ba as our heliocentric reference frame, then added Bb by specifying its orbital parameters with respect to Ba and shifted the reference system to the centre of mass of BaBb. Later, we used the AB orbital parameters to simulate a third body with a combined mass M_Aa + M_Ab which corresponds to the centre of mass of the A system. We then saved the initial 3D positions x₀ and velocities v₀ of this third body ‘AaAb’ using the centre of mass of BaBb as the reference frame. We then set up a new simulation, only for the AaAb system to get the initial positions of Aa and Ab, x_Aa, 0, x_Ab, 0 and velocities v_Aa, 0, v_Ab, 0 with respect to the centre of mass of the AaAb pair. All the positions and velocities for all four stars were calculated at the same reference time, in our case we used T₀ of the AB orbit. Finally, we set up the final simulation by adding Ba, followed by Bb by specifying its orbital parameters with respect to Ba, which moved to the centre of mass of BaBb. We then added Aa by specifying its initial position and velocity calculated earlier, the position and velocity are x₀ + x_Aa, 0 and v₀ + v_Aa, 0, respectively, and we then did the same for Ab.

Figure C.1 shows the positions of the four stars as we integrated the simulation in time for both solutions I and II, overlapped with the location of the disk. The centre of mass of BaBb is located at (0,0).

Fig. C.1. Integrated orbits at the times of transit of AaAb behind the disk surrounding BaBb, using the best-fit parameters. The disk and the four orbits are referred to the centre of mass of BaBb located at (0,0).

Appendix D: Flux ratio estimation

We used evolutionary track from Baraffe et al. (2015), assuming an age of 10 Myr, and synthetic photometry with a BT-Settl model grid, provided by the Spanish Virtual Observatory (SVO) web service¹³ to estimate the flux ratio corresponding to the dynamical masses obtained in this work. The theoretical flux from the BT-Settl model was scaled by the multiplicative dilution factor M_d = (R/D)², R being the stellar radius and D the distance to the observer (see Tables D.1 and D.2).

All Tables

All Figures

	Fig. 1. Schematic view of HD 98800 orbital configuration. The BaBb subsystem hosts a circumbinary disk in polar configuration and is orbiting around the AaAb binary in a highly inclined orbit with a semi-major axis of ≈45 au.
In the text

	Fig. 2. Squared visibility and closure phase measurements from one observation of AaAb taken in March 2021. The data are in blue, while the red dots represent the best binary model fitted with CANDID for this epoch. Bottom panels: residuals in the number of sigmas.
In the text

	Fig. 3. Detection level map from CANDID for the observation of AaAb taken in April 2019. The colourbar shows the significance of the companion detection in the number of sigmas. The red cross points to the best-fit position.
In the text

	Fig. 4. Diagram of the orbit of the secondary star around the centre of mass (yellow plane) and the reference plane (grey). This diagram follows the orbital convention of exoplanet.
In the text

	Fig. 5. Phase-folded RVs orbit for BaBb. The systemic velocity γ for each set of observations was subtracted. The solid line corresponds to the best-fit model. Upper panel: RVs of Ba, and lower panel: Bb.
In the text

	Fig. 6. Best orbital solution for BaBb. The solid line corresponds to the best-fit model and the shaded area to the 1σ region. The primary star Ba is located at the origin. The relative positions of Bb are plotted as filled dots. The error ellipses from PIONIER astrometry are smaller than the markers.
In the text

	Fig. 7. Posterior distribution of masses and parallax of AaAb subsystem assuming the HIPPARCOS value (red) and the solution from BaBb fitting (blue) as a prior distribution of the parallax in our MCMC model. The red and blue lines highlight the median of each distribution.
In the text

Fig. 8. Best orbital solution for AaAb. In both panels, the solid line corresponds to the best-fit model. Bottom panel: the primary star Aa is located at the origin. The relative positions of Ab are plotted as filled dots; the error ellipses from PIONIER astrometry are smaller than the marker. Upper panel: the coloured markers correspond to the primary star RV measurements. The systematic velocity γ for each set of observations was subtracted.

In the text

	Fig. 9. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. The solid black line corresponds to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue one assumes large uncertainties (σ ∼ 0.1″). The red dots correspond to the astrometric measurements. For better visualisation, the error bars are shown only in Fig. 10.
In the text

Fig. 10. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. In all panels the solid line corresponds to the best-fit model and the shaded area to the 1σ region. The solid black lines correspond to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue ones assume large uncertainties (σ ∼ 0.1″). The red dots correspond to the astrometric measurements. The error bars shown in the astrometry before 1991 correspond to the large uncertainty case.

In the text

Fig. 11. Best orbital solution for AB outer orbit for both uncertainty assumptions in the astrometry before 1991. In both panels the solid line corresponds to the best fit model and the shaded area to the 1σ region. The solid black lines correspond to the best orbital solution assuming small uncertainties (σ ∼ 0.02″) and the blue ones assume large uncertainties (σ ∼ 0.1″). The dots markers correspond to the RV measurement of systemic velocities from our orbital solutions and the one obtained from CO modelling.

In the text

	Fig. 12. Mean MEGNO value for the 12 × 10 matrix for different integration times. Error bars correspond to 1σ.
In the text

Fig. 13. Probability density plot of 1000 realizations of the light curve for the occultation of AaAb behind the disk surrounding BaBb for solutions I and II (top and bottom, respectively). The colour bar shows the probability of getting a determined flux at a given time, such that the sum along each of the columns is normalised to unity. In orange we show the light curve for the best-fit parameters (Table 4).

In the text

	Fig. B.1. Squared visibilities from Keck Interferometer observations published in Boden et al. (2005). The black circles represent the observed values and the red crosses represent the best-fit BaBb binary model from this work.
In the text

	Fig. B.2. Posterior samples of AaAb orbital parameters. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.
In the text

	Fig. B.3. Posterior samples of BaBb orbital parameters. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.
In the text

	Fig. B.4. Posterior samples of AB orbital parameters for solution I. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.
In the text

	Fig. B.5. Posterior samples of AB orbital parameters for solution II. Contoured sub-panels show the distribution of points from the MCMC chains, where high-density regions are indicated by the greyscale and contours. Histogram sub-panels show the posterior distributions, with median and 68% confidence intervals marked by dashed lines, with titles quantifying those ranges.
In the text

	Fig. C.1. Integrated orbits at the times of transit of AaAb behind the disk surrounding BaBb, using the best-fit parameters. The disk and the four orbits are referred to the centre of mass of BaBb located at (0,0).
In the text