© S. A. Zarrilli et al. 2022

1. Introduction

The masses and ages of young stellar objects (YSOs) are commonly derived from a comparison of observed positions on a colour-magnitude diagram (CMD) with theoretical models (e.g., Siess et al. 2000; Baraffe et al. 2015; Choi et al. 2016). Evolutionary tracks of YSOs are very sensitive to mass and as such need to be calibrated from observed systems with well-constrained masses. The relative paucity of higher-mass Herbig Be YSOs (≳5 M_⊙) compared with their lower-mass T Tauri counterparts means that models of higher-mass stars are less thoroughly calibrated. The high effective temperature of Herbig Be stars (≳15 000 K), as well as their often-uncertain ages due to mass loss from stellar winds, makes it much less straightforward to calculate their masses (Massey et al. 2012). Stassun et al. (2014) find that predicted and measured masses can differ by ∼10%.

The gold standard for deriving model-independent masses is to take advantage of the orbital mechanics of binary systems. If both astrometric and radial velocity (RV) data are used, the derived orbital parameters can be combined with reliable distance estimates (if any exist) to extract dynamical masses for the individual objects. This requires observation of a binary system at multiple epochs spread over a substantial fraction of the orbit, so targets with relatively short orbits and small separations are the best candidates. The need for precise astrometry on very small angular scales (∼1 mas) has hugely benefitted from the relatively recent development of optical and near-infrared (NIR) interferometry, which has ‘unlocked’ a larger tranche of suitable systems compared to even 15 years ago.

A typical feature of YSOs is the presence of substantial amounts of circumstellar material left over from stellar formation, which takes the form of a disc due to conservation of angular momentum. In single systems, dispersal of the disc occurs from a combination of accretion onto the star, depletion from stellar wind, and condensation into protoplanets, with a typical disc lifetime of 1–3 Myr (Li & Xiao 2016). However, the picture is more complicated when binary systems are concerned. If the binary is widely separated, each individual star can host its own circumstellar disc, but for many luminous and close Herbig Ae and Be binaries, a single circumbinary disc is the only possible structure for circumstellar material (Pichardo et al. 2005). This is due to dynamical interactions between the stars and the disc, which can affect both the accretion properties of the system and the disc’s shape and lifetime, although it is unclear under what conditions they will either delay or accelerate disc dispersal (Cieza et al. 2009). Dynamical truncation will affect the potential of the disc to form planetary systems by removing or rearranging the material available for planet formation. The known population of circumbinary planets has been substantially increased by Kepler observations (e.g., Doyle et al. 2011), and numerical simulations suggest that features rare in planets around single stars, such as large eccentricities and planet-star misalignment, are more common in circumbinary systems (Chen et al. 2019). Further studies on well-characterised young multiple star systems are also essential for studying other dynamical mechanisms that might shape the architecture of exoplanetary systems, for instance by moving disc material onto oblique orbits (Kraus et al. 2020).

The focus of this study is the multiple system MWC 166 (also known as HD 53367, HIP 34116, and V750 Mon). This is a hierarchical triple system, with a close spectroscopic binary (MWC 166 A) orbited by a wide companion (MWC 166 B) at a separation of 0.6″ (Fabricius et al. 2002). The RV variations of the spectroscopic binary were first reported by Finkenzeller & Mundt (1984), and tentative spectroscopic orbit solutions have been presented by Corporon & Lagrange (1999) and Pogodin et al. (2006). In this paper we focus on this inner spectroscopic binary, the components of which have been labelled MWC 166 Aa and MWC 166 Ab throughout.

MWC 166 is located in the nearby OB association Canis Major OB1, whose age has been estimated to be ∼3 Myr (Clariá 1974). The object also features significant mid-infrared to millimetre excess, which is indicative of a disc around MWC 166 A. The distance to the OB association has been estimated to be 1150 ± 140 pc (Clariá 1974). Subsequent photometric measurements taken from a larger number of sources broadly agree with this value and have more tightly constrained it to 990 ± 50 pc (Shevchenko et al. 1999; Kaltcheva & Hilditch 2000).

Here, we present NIR interferometric observations obtained with the Very Large Telescope Interferometer (VLTI) and the Center for High Angular Resolution Astronomy (CHARA; ten Brummelaar et al. 2005) Array, which have allowed us to derive a first astrometric orbital solution of the system. We present our observations in Sect. 2, which is followed by a discussion of our modelling approach (Sect. 3). We derive the orbit solution and dynamical mass constraints in Sect. 4, and spectral line analysis results are presented in Sect. 5. A discussion on the distribution of circumstellar material – both in the dust continuum and in prominent K band emission lines – is presented in Sect. 6, and our conclusions are summarised in Sect. 7.

2. Observations

Near-infrared interferometric observations were taken over a period of 8 years, mainly using the PIONIER (Le Bouquin et al. 2011) and GRAVITY (GRAVITY Collaboration 2017) four-beam combiners at the VLTI. All VLTI observations employed the 1.8-metre Auxiliary Telescopes. Longer baselines were provided by four-telescope observations using the CHARA Array instrument MIRC-X (Kraus et al. 2018; Anugu et al. 2020).

The GRAVITY observations were taken in the K band (1.99–2.45 μm) as part of ESO programme 098.C-0910(A). GRAVITY observations include data from the fringe tracker, which operates at a spectral resolution of ℛ = Δλ/λ ∼ 22, as well as the science combiner either in medium (ℛ ∼ 500) or high (ℛ ∼ 4000) resolution (GRAVITY Collaboration 2017). Our observations achieved an angular resolution up to (λ/2B_max) = 1.6 milliarcseconds (mas) on the longest baselines (B_max = 130 m in length), sufficient to spatially resolve the components of MWC 166 A (GRAVITY Collaboration 2017). The reduction pipeline used was the GRAVITY data reduction pipeline¹ running in the ESOreflex v2.9.1 environment (Freudling et al. 2013). Besides the statistical uncertainties computed by the GRAVITY pipeline, we include 5% and 1° errors, for the visibility and closure phase, respectively, to account for calibration uncertainties.

The PIONIER observations covered the H band (1.59–1.75 μm) and were obtained as part of multiple ESO programmes: 102.C-0701, 104.C-0737, and 106.21JU. These data were recorded over 6 channels at spectral resolution ℛ ∼ 40. The reduction pipeline used was pndrs v3.52 (Le Bouquin et al. 2011). We also included published data from the large programme 190.C-0963² (e.g., Lazareff et al. 2017; Kluska et al. 2016), over 3 spectral channels and with a resolution of ℛ ∼ 15.

MIRC-X was used in its H band mode as part of the programme 2020B-M7, using four of the six CHARA Array telescopes. CHARA’s much longer maximum baseline of 330 m allowed us to probe the object geometry at a nearly three-times higher resolution than possible with VLTI. The MIRC-X v0.9.5 pipeline³ (Anugu et al. 2020, Sect. 4) was used to reduce the data.

A full description of the observations is provided in Table 1. Each observation was calibrated by observing suitable calibrator stars with known uniform disc diameters (taken from Bourges et al. 2017), to account for atmospheric absorption and instrument response. The calibrators were also inspected for signatures of binarity to ensure only single stars were observed. In order to improve the (u, v)-coverage, we grouped the individual measurements into epochs. However, this proved to be difficult due to the rapidly changing orbit of the system. Based on the literature RV orbital period of ∼183 days (Pogodin et al. 2006), a variation of about one degree in position angle (PA) per day is to be expected. Considering that the PA uncertainties in our binary model fits are on the same order of magnitude as this (see Table 3), each epoch should include data from at most two consecutive nights, ensuring that the relative positions of the two components do not change significantly during each epoch. The observations for which this consolidation was performed are marked accordingly in Table 1.

3. Modelling

3.1. Continuum modelling of the system

We fit the interferometric visibility and closure phase data at each epoch using the Exeter in-house geometric modelling pipeline (e.g., Kreplin et al. 2018). As the stellar radii of MWC 166 Aa and Ab are expected to be ∼0.04 mas at the distance calculated by Kaltcheva & Hilditch (2000), we assume that the stellar photospheres can be modelled as point sources.

Initially, the visibilities and closure phases of MWC 166 A were fitted with the following free parameters: separation (ρ); PA⁴ of the secondary component from the primary (θ); and the flux contribution of the secondary to the total flux in the model (f₂/f_tot). The primary flux contribution was kept fixed

Due to the changes in the f₂/f_tot flux ratio, the secondary is at some epochs brighter than the primary in our NIR wavelength bands. While there is some variability to the system as a whole over year-length timescales (Pogodin et al. 2006), the relative brightness of the two stars has not been previously recorded, so this was an unexpected finding. In light of this, we restricted θ either to the [0°, 180°] or [180°, 360°] range, where the quadrant was chosen for each epoch to achieve an astrometric orbit that is consistent with the spectroscopic orbit of Pogodin et al. (2006). This was done to ensure that the primary and secondary components of the system were correctly identified at each epoch.

Initially, we did not account for contributions from possible dust emission, consistent with the low measured infrared excess emission in the K band (Tjin A Djie et al. 2001). The two-point-source model fits allowed us to derive the astrometry of the two components of the system, but yield a flux ratio that changes significantly between epochs, both in the H and K band. These models also consistently over-predicted the visibilities, as can be seen from the red points in Fig. 1, leading to large reduced χ² values, in particular on the visibilities (e.g., ${\chi }_{\text{vis}}^2$ > 16 for the MIRC-X data). In order to reduce this systematic error, we also conducted fits that include extended emission.

Fig. 1. Visibilities (and associated residuals) of MIRC-X models. At V ≳ 0.6, the purely point-source model overshoots the observed data points substantially.

3.1.1. Evidence for extended circumbinary disc emission

We tried modelling the extended flux assuming three different geometries: a Gaussian with full width at half maximum (FWHM) σ, seen under inclination i and a major axis (east of north) PA Θ_ext; a ring with radius R and a thickness of 0.2R, seen under inclination i and a PA Θ_ext; and as over-resolved flux ‘background’ (modelled as a circular Gaussian with σ = 1000 mas, i.e. filling the field-of-view). The integrated flux contribution of the extended emission component to the total flux in the model is f_ext/f_tot, where we define f_tot ≡ f₁ + f₂ + f_ext.

For the VLTI epochs, the different geometries for the extended emission returned improved χ² values over the two-point-source model, but no one extended model had consistently smaller χ² values over all epochs. Adopting a Gaussian geometry for the extended emission component results in ${\chi }_{\text{vis}}^2$ values ranging from 0.34 to 2.24, while adapting a ring geometry results in ${\chi }_{\text{vis}}^2$ = 0.49...4.65, and over-resolved flux results in ${\chi }_{\text{vis}}^2$ = 0.57...6.81. For comparison, the pure point-source model has ${\chi }_{\text{vis}}^2$ between 0.94 and 31.86. Closure phase χ² values were found to be almost completely model-independent, with the models returning ${\chi }_{\text{CP}}^2$ <  2.84 (point-sources), ${\chi }_{\text{CP}}^2$ <  1.62 (background), ${\chi }_{\text{CP}}^2$ <  1.13 (Gaussian), ${\chi }_{\text{CP}}^2$ <  1.43 (ring).

The MIRC-X data probes roughly three-times higher spatial frequencies than the VLTI data. Data taken over a larger range of spatial frequencies allows us to probe farther lobes of the visibility curve, which helps to more reliably distinguish the effects of the extended emission from the sinusoidal binary modulation. Figures 1 and 2 respectively show the resultant model visibilities and closure phases obtained from our geometric models described above, overlaid over the data. While the closure phases are well described by the original two-point-source model independently of any extended flux or lack thereof, for the visibilities this is not the case. The two-point-source model, represented by the red points in Fig. 1, clearly overpredicts the visibility compared to the other models, especially in the high-visibility regime (where V ≳ 0.6). The model parameters corresponding to Figs. 1 and 2 are shown in Table 2.

Fig. 2. Closure phases (and associated residuals) of MIRC-X models.

By examining the χ² values both for visibility and closure phase in Table 2, it can be seen that the background model provides a significant improvement on the two-point-source model. A ring profile provides a similar, or even slightly better fit (see Table 2), but introduces three additional free parameters while providing only a marginal improvement in the goodness of the fit. We also conducted a fit for a Gaussian model (σ, i, Θ_ext), which returned similar χ² values as the ring model, but we found that the parameters i and Θ_ext did not converge to a value independent of the boundary conditions chosen, while the flux parameters f₂/f_tot and f_ext/f_tot were found to be consistent with the background model. As such, we favour the background model over the Gaussian model.

Therefore, we adopt the over-resolved background model as geometry for the extended emission for all epochs and instruments, likely representing scattered light from the disc. This minimised the model complexity and degrees of freedom. The relative astrometry was therefore fitted with four free parameters: separation (ρ); PA of the secondary component from the primary (θ); secondary flux as fraction of the total flux (f₂/f_tot) and extended flux as fraction of the total flux (f_ext/f_tot).

It is apparent that the relative astrometry of the binary (ρ, θ) does not depend much on whether extended flux is included in the fit. The value of f₂/f_tot and f_ext/f_tot change depended on whether extended emission is included in the fit, but is rather independent of the geometry of the emission.

3.2. Modelling of the K band He I and Br γ lines

For the GRAVITY data, we simultaneously recorded high-resolution data (Δλ/λ = 4000) for all epochs, in addition to the low-resolution data used to establish the relative astrometry. As can be seen from MWC 166 A’s K band spectrum, there are strong line features at 2.058 and 2.166 μm (Fig. 3), corresponding to helium-I and Brackett-γ emission, respectively.

Fig. 3. Continuum-normalised spectra around the He I and Brγ lines. The different epochs have been offset for clarity. Epoch-dependent variations are visible. Dates are given in the format YYYY-MM-DD.

In order to model the spectral lines, we used the fitting tool PMOIRED⁵, which is capable of fitting closure phases, differential phases, differential visibilities, and line spectra simultaneously, both over the continuum and over specific spectral windows. After finding the continuum geometry of the system for each epoch (Sect. 3.1.1), we introduced new model components to fit the He I and Br γ lines individually, using geometries of varying complexity. We initially used a single Gaussian of FWHM σ = 0.1 mas and left the position of the line-emitting region as a free parameter. At all epochs, this resulted in only very small spatial displacements from the origin, suggesting the primary component is responsible for the majority of the emission in the system. However, this singular emission zone near the primary is perhaps too simplistic a model. If we examine the flux intensity of the two lines (Fig. 3), there are indications of time-dependent variability, as well as signs of double-peaked lines at several epochs. This could imply that the emission originates from both stars, or it could a signature of the gas kinematics. Furthermore, the differential phases over the spectral lines show signatures of rotating emission that are not accurately modelled by the single-Gaussian model. These features would be most naturally explained by a circumprimary disc.

3.2.1. Circumprimary gas disc model

The circumprimary disc model we used is based on the model described in Frost et al. (2022), which was used to model the binary system HR 6819. The model comprises several components and has 12 free parameters that describe the entirety of the system.

Firstly, the stars themselves are modelled as uniform discs with diameters corresponding to twice the stellar radii we obtained from our continuum fit (see Sect. 4.4), and the secondary component is given a displacement from the primary’s position at the origin, as well as a continuum flux f₂/f₁, while the primary flux is fixed to f₁ ≡ 1.

The line emission is subsequently modelled to be originating from two regions, one blueshifted and the other redshifted, to represent the approaching and receding part of a rotating disc (labelled B and R, respectively). Each of these components was given its own spatial displacement (x_i, y_i), with the size of the emitting region following a Gaussian profile with FWHM ${\sigma }_i=\frac{1}{2}\sqrt{x_i^2+{\mathit{y}}_i^2}$. We additionally modelled the line components in the spectral domain. Each component is given a flux profile F_i = f_i + F_L, consisting of a Lorentzian component F_L, which was kept equal for both wings, and a flat component f_i, accounting for the differences in line strength between the components (which can be seen to vary by epoch in Fig. 3). Each line wing is centred on a wavelength that is displaced from the central line wavelength λ₀, such that λ_B = (λ₀ − Δλ) and λ_R = (λ₀ + Δλ).

The resultant fitted parameters of the model described above are presented for all GRAVITY epochs in Sect. 5. The model is fitted simultaneously to the telluric-corrected spectrum, closure phase, differential phase, and differential visibility.

4. Results: Orbital solution and mass–distance constraints

4.1. Binary astrometry and orbital fit

The parameters of our best-fit continuum model with background component are listed in Table 3. The model visibilities that correspond to the best-fit model are shown in Fig. 4, while closure phases are shown in Fig. 5. Our modelling script makes use of Markov chain Monte Carlo module emcee (Foreman-Mackey 2016) to explore the parameter space and obtain error estimates from the posterior probability distribution. We show the corner plot for a representative epoch (2 February 2018) in Fig. 6. An important source of systematic uncertainty that affects primarily the derived separations is the wavelength calibration, and we account for this by including a systematic uncertainty of 5% for the separations listed in Table 3.

Fig. 4. Observed continuum visibilities (black) and corresponding models (red) plotted against spatial frequency for all epochs. The model fit includes extended background emission. Dates are given in the format YYYY-MM-DD.

Fig. 5. Observed continuum closure phases (black) and corresponding models (red) plotted against spatial frequency for all epochs. The model fit includes extended background emission. The scaling on the vertical axis was adjusted for each epoch. Dates are given in the format YYYY-MM-DD.

Fig. 6. Corner plot showing the possible correlations between free parameters (ρ, θ, F2 = 100 ⋅ f2/ftot and F3 = 100 ⋅ fext/ftot, respectively) for epoch 6 February 2018 of the continuum GRAVITY data, where the model fit includes extended background emission.

Using the relative astrometry for each epoch, we fitted a Keplerian orbit using the standard Campbell elements: P = orbital period; T₀ = epoch of periastron passage; a₁ = semi-major axis of primary component; i= orbital inclination (to the line of sight); e = eccentricity; Ω = longitude of ascending node; ω= longitude of periastron; K₁ = orbital curve semi-amplitude of primary component; V₀ = RV of the system’s centre of mass.

We used two fitting approaches:

ORBITX code. This code⁶ (Tokovinin 1992) fits orbits using both astrometric and RV data simultaneously. We modified the ORBITX code to account for uncertainties in both ρ and θ, a feature absent from the original code, which only uses uncertainties on ρ.

Grid-search algorithm. We used the grid-search algorithm developed by Kraus et al. (2009) to construct a grid of orbital solutions in the P, T₀, and e parameter space, where the remaining elements a, i, Ω, ω are determined from the Thiele-Innes elements. We explored the parameter space around P = 0.480...1.100 yr, T₀ = 2019.5...2020.5 yr (step sizes of 0.001 yr), and e = 0.100...0.600 with a step size of 0.001, and selected the solution with the lowest combined residuals in RV and astrometry. We then repeated the process with smaller step-sizes around the initial solution (a factor of ten for all parameters, = 0.0001 yr, 0.0001 yr, 0.0001, respectively, for P, T, e), in order to increase precision. Uncertainties were calculated by examining the χ² curve for each parameter.

The best-fit orbit solutions found with these methods are listed in Cols. (3), (4) of Table 4, and are overplotted on the data in Figs. 7 and 8. Both the orbits provide a very good fit to the existing data, and provide similar results for all parameters. This is despite a substantial portion of the orbit still lacking astrometric observations. We adopt the orbit from Col. (3) when discussing derived quantities in the subsequent sections.

Fig. 7. Astrometric orbit solutions derived using both the ORBITX code (blue line) and the grid-search code (red line). The primary star is kept fixed at the origin, and the x and y axes show displacement in right ascension and declination, respectively. The dotted lines connect the ascending and descending nodes of each orbit.

Fig. 8. Radial velocity measurements of MWC 166 A taken in the period 1994–2005 plotted against orbital phase. The blue and red fitted curves correspond to the orbits specified in Cols. (3) and (4) of Table 4, respectively. The dotted black line shows the velocity of the system’s centre of mass (V0) for the grid-search method, and the solid grey line shows V0 for the ORBITX orbit.

4.2. Comparison to RV orbit

Corporon & Lagrange (1999) and Pogodin et al. (2006) derived orbits for MWC 166 A from the RV data, with the more recent of the two being a refinement including additional RV points. The orbital parameters for this RV orbit are shown in Col. (2) of Table 4. Our spectroscopic and astrometric orbital solutions differ substantially from the earlier RV-only orbit. The most notable difference is in the orbital period, which we calculated as almost exactly twice the length of the Pogodin et al. (2006) orbit. This doubling of the period was only discernible thanks to our astrometric data, as using the RVs alone provides an equally good fit to both orbits. We also found the orbit to be much more elliptical than previously thought, with its eccentricity of 0.498 ± 0.003 being much larger than that of the RV orbit (e =  0.28 ± 0.03). A newly determined parameter from our orbit is the inclination, with a value of i =  53.6 ± 0.3°.

4.3. Dynamical system mass

According to Kepler’s third law, the period of an orbit P is proportional to the cube of the semi-major axis of the orbit a. Using the usual angular diameter-distance relation a [au] = a [″] × d [pc], we can show the dependence of the total system mass M_tot ≡ (M₁ + M₂) (in solar masses) on distance d (in parsecs), where G is the gravitational constant and M_1, 2 the masses of the two stars:

$$\begin{array}{c}\hfill M_{\mathrm{tot}}=\frac{4\pi a^3\mathcal{Z}}{G{P}^2}{d}^3,\end{array}$$(1)

where $\mathcal{Z}=1684.14{\text{m}}^3{M}_{\odot }{}^{-1}$ is a constant introduced to account for the change in units from au to metre, and from kg to M_⊙. In the above equation, the total mass is also known as the ‘dynamical mass’, signifying it is derived from fitting the dynamical orbit of the system. Since M_tot ∝ d³, a reliable distance value is needed to obtain rigorous mass estimates.

Unfortunately, there are several conflicting distance estimates for this system in the literature. As mentioned in Sect. 1, photometrically calculated distances have placed MWC 166 at a distance of roughly one kiloparsec. Parallax observations have corroborated the photometric distances for a majority of the other individual members of CMa OB1, but are not consistent in the case of MWC 166. HIPPARCOS (ESA 1997) measured a distance of 247 ± 82 pc, while Gaia Data Release 2 (DR2) parallaxes correspond to an even shorter distance of ${131}_{-13}^{+16}$ pc (Bailer-Jones et al. 2018) – roughly ten times closer than the ∼1 kpc to its parent association. This discrepancy can likely be explained by the binarity of the system, as DR2 does not solve for source multiplicity. Indeed, the recent release of preliminary results from Gaia Early Data Release 3 has brought the parallax distance to MWC 166 closer to the photometric distance values, albeit with very large uncertainties (d ∼ 1600 ± 700 pc, Bailer-Jones et al. 2021). Our Keplerian mass-distance relation (Eq. (1)) shows that distances of ≲650 pc correspond to masses of M_tot < 5 M_⊙, clearly below the mass threshold for MWC 166 Aa’s spectral type of B0III (Fairlamb et al. 2015). Conversely, the photometric distances offer more physically realistic values – the most recent distance estimate of 990 ± 50 pc returns a system mass of M_tot = 17.1 ± 2.7 M_⊙, which is also in agreement with the prediction of 20–25 M_⊙ by Pogodin et al. (2006). Additionally, any parallax distances are unreliable due to having been calculated with the assumption of a six-month orbit, which we have shown is too short by a factor of two. In light of these points, in this work we have treated the photometric distances preferentially. We have therefore taken the literature distance to MWC 166 to be the most recent photometric distance, 990 ± 50 pc (Kaltcheva & Hilditch 2000).

4.4. From combined mass to individual masses and other properties

Our full orbital solution allows us to constrain model-independent individual masses for the first time. If the system’s orbital period (P), eccentricity (e) and inclination (i) are known, as well as the RV semi-amplitude of the primary component (K₁), it is possible to calculate the binary mass function f (e.g., Boffin 2012; Curé et al. 2015):

$$\begin{array}{c}\hfill f\equiv \frac{K_1^{3}P(1-e^2{)}^{3/2}}{2\pi G}=\frac{{(M_{2}sini)}^3}{M_{\mathrm{tot}}^2},\end{array}$$(2)

which equates the mass of the secondary component to the other elements. Rearranging for M₂ gives

$$\begin{array}{c}\hfill M_2=\frac{1}{sini}[\frac{K_1^{3}P(1-e^2{)}^{3/2}}{2\pi G}·M_{\mathrm{tot}}^2{]}^{1/3},\end{array}$$(3)

and the mass of the primary can therefore be trivially found through M₁ = M_tot − M₂.

Using this method we determined the masses of MWC 166 Aa and MWC 166 Ab as M₁ = (12.19 ± 2.18) M_⊙ and M₂ = (4.90  ±  0.52) M_⊙, respectively. Our M_dyn measurements for both components were used to derive the remaining stellar parameters, and their respective confidence bounds, from theoretical evolution tracks. We used CMD 3.6⁷ to generate a solar metallicity (Z = 0.0152) isochrone table from PARSEC 1.2S (Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014; Marigo et al. 2017; Pastorelli et al. 2019).

The evolutionary status of MWC 166 A has not been conclusively established. Analysis of the spatial distribution of O- and B-type stars in its parent OB association, CMa OB1-R1, resulted in an estimate for its age of ∼3 × 10⁶ yr (Clariá 1974). At this age, MWC 166 Aa is likely already onto the main sequence, while MWC 166 Ab might still be in its pre-main-sequence stage (Tjin A Djie et al. 2001). Later work by Herbst & Assousa (1977) showed that CMa OB1-R1’s main stars are located on the rim of an expanding shell of neutral hydrogen, consistent with star formation being triggered by a supernova within the last ∼500 kyr. Due to the uncertain age of the system, we selected these two potential ages, as well as a reasonable intermediate age (1 × 10⁶ yr), upper bound (1 × 10⁷ yr), and lower bound (1 × 10⁵ yr). Then, for each of the estimated ages, we drew an isochrone of an appropriate age from the CMD table and performed a quadratic interpolation over the mass points within it to find the value for each parameter that corresponds to M_dyn at the given age. This process was repeated for the upper and lower bounds on M_dyn to find the confidence bounds for each parameter at the given age. Our estimates for the stellar parameters using this method are presented in Table 5. Figure 9 shows the isochrones for each age estimate plotted on a Hertzsprung-Russell diagram of log(L) versus log(T_eff), as well as the interpolated values for the primary and secondary components as large circles and triangles, respectively.

Fig. 9. Hertzsprung-Russell diagram showing PARSEC 1.2S isochrone tracks for each age. Superimposed are the interpolated L, Teff values for the primary star (large circles) and secondary star (triangles), for all ages. The black cross shows the location of the primary’s parameters as determined by Fairlamb et al. (2015).

From Fig. 9, it can be seen that, at the ages of 500 kyr, 1 Myr, and 3 Myr, the primary star is on the main sequence, with negligible variations in parameters between these ages. At the upper bound age of 10 Myr, the primary is in the process of beginning to evolve beyond the main sequence. Conversely, it is clear from the Hertzsprung-Russell diagram that the secondary component is predicted to be on the main sequence for all ages save the youngest, where it is still in its protostellar stage. Also in Fig. 9 is shown the location of the primary as calculated by Fairlamb et al. (2015), which appears to agree with all ages except the upper bound of 10 Myr. In Sect. 6.4, we generalise this process to all generated isochrones to attempt to constrain the age of the system.

5. Results: Modelling the gas distribution and kinematics in the He I and Br γ lines

The circumprimary disc model described in Sect. 3.2 produces an excellent fit to the data at all epochs.

Figure 10 shows the results of the fit for the epoch 14 March 2017, where panel a shows the spectra, closure phases, differential phases, and visibilities measured with GRAVITY overplotted with prediction from the best-fit model. The three left-most plots in panel b show synthetic images computed from the best-fit model for three representative wavelengths. The central panel shows the brightness distribution at λ₀, while to the left and right of this the wavelength is redshifted and blueshifted respectively by a factor of λ₀ ± 3Δλ. The right-most panel of the figure shows the flux contribution of each synthetic model component.

Fig. 10. Results of the Br γ line modelling, for the epoch 14 March 2017. (a): (u, v)-coverage for the observations associated with the epoch, coloured by baseline pair. (b): telluric-corrected normalised flux (labelled NFLUX, black lines), overplotted with flux in the best-fit model (red line). (c)–(e): data from each GRAVITY exposure (black lines), overplotted with quantities computed from best-fit model (red lines). The observables are closure phase for each telescope triplet (T3PHI), differential phase for each baseline (DPHI), and visibility for each baseline (|V|), respectively. (f)–(h): brightness distribution corresponding to the best fit, for three representative wavelengths. (i): synthetic line strengths and profiles for the two spectral components in the model (red and blue), in addition to the monochromatic continuum flux associated with the primary component (silver) and the secondary component (gold). The parameters corresponding to the fits can be found in Table 6.

The line modelling results for all epochs (analogously to Fig. 10) are shown in the Appendix, and the findings are summarised in Table 6, for both spectral lines of interest.

The results show substantial evidence of line emission being localised around the primary star, consistent throughout all four GRAVITY epochs. The redshifted component (R) is displaced consistently towards the north-west of the primary star, while the blueshifted component (B) is located south-east of the primary. The displacement of the two components from the star is generally (∼0.25 mas). The average displacement for He I was found to be (10.5 ± 0.4)R₁, and for the Br γ this was (11.5 ± 0.4)R₁. The derived values for each epoch can be seen in Col. (5) of Table 7. The relative intensities of the B and R components also seem to be variable by epoch and by spectral line. From Cols. (3) and (4) of Table 7, it can be seen that, for the He I line, the red wing contributes 60% of the flux for the first two epochs, a similar level of intensity as the blue wing during the latter two epochs. On the other hand, the Br γ line flux tends to be more equally distributed between the two components, apart from in the epoch 6 February 2018, where the blue component is substantially stronger. The emission of the redshifted and blueshifted line wing are displaced along an average PA of 134.0 ± 1.1°, indicating that this is the sky-projected PA of the major axis of the rotating disc. Accordingly, the rotation axis of the discs, and likely also the primary star, is oriented along PA 44.0 ± 1.1° on sky.

6. Discussion

6.1. Evidence for circumbinary dust

Our interferometric observations in the continuum allow us to quantify the excess emission through our visibility fit (Sect. 3.1.1), where we find that ∼2% and ∼5% of the excess emission are associated with extended flux in the H band and K band, respectively. The H band observations also show ‘spikes’ in the extended flux contribution in the epochs 24 December 2019 and 19 November 2020, with f_ext/f_tot contributing up to 15% of the total flux for these epochs. The geometry of the extended dust is still poorly constrained, but our modelling shows the emission originates from scales at least twice larger than the binary separation vector (for a ring model), but likely even larger (for a Gaussian or background flux model; Table 2 and Sect. 3.1.1). This is comparable to the dynamical truncation radius predicted by Artymowicz & Lubow (1994) for circular binaries (∼1.7a).

The expected Silicate dust sublimation radius (R_s) can be estimated with

$$\begin{array}{c}\hfill R_{\mathrm{s}}=1.1\sqrt{Q_{\mathrm{R}}}(\frac{L}{1000L_{\odot }}{)}^{1/2}(\frac{T_{\mathrm{s}}}{1500\mathrm{K}}{)}^{-2}\mathrm{au},\end{array}$$(4)

where L is the bolometric luminosity of the star(s) irradiating the disc and T_s is the dust sublimation temperature (Monnier & Millan-Gabet 2002). Q_R ≡ Q_abs(T_*)/Q_abs(T_s) is the ratio of dust absorption efficiencies for radiation for the incident and re-emitted field, which we fix to Q_R = 1 in order to estimate the rim location for large micron-sized dust grains. The above equation has several simplifications compared to the physical system, chief among which is the assumption of a single star rather than a binary. While dust sublimation radii have been calculated for binary systems (e.g., Nagel et al. 2010), dynamical interactions cause the inner dust rim to be at a larger radius than R_s. A detailed calculation is therefore outside the scope of this work, but Eq. (4) can still give an order-of-magnitude estimate for the minimum dust inner rim radius of the circumbinary disc, assuming a luminosity equivalent to that of the two components combined. Substituting the sum of the luminosity values from the central age estimate of 1 Myr in Table 5 gives R_s = 2.6 au, 4.6 au and 7.1 au for dust sublimation temperatures of 2000 K, 1500 K and 1200 K, respectively. This suggests that individual dusty circumstellar discs are likely not present in the MWC 166 A system. Due to the calculated semi-major axis of the orbit being comparable to the smallest of the above estimates (a₁ = 2.61 au at d = 990 pc), it is expected that substantial individual dusty discs around each of the components are not capable of surviving for extended periods of time due to dynamical interactions (e.g., Mathieu 1994). Accordingly, we associate the extended continuum emission with a circumbinary disc component instead of circumstellar disc component.

6.2. Evidence for variable extinction or circumstellar material

The relative flux contribution of the two point sources in our model was found to be variable, especially in the H band continuum. The flux associated with the secondary appears to increase around phase ∼0.9, and was found to be as bright as the primary at one epoch (24 December 2019) and even brighter than the primary at two epochs (16 December 2019 and 28 December 2020). These last two epochs were both taken on the compact AT configuration, and returned visibilities very close to unity – as well as closure phases close to zero – at all probed baselines (see Figs. 4 and 5). This, and the more limited (u, v)-coverage caused by the shorter baselines, can result in ambiguities in the model (Anthonioz et al. 2015). To rule this out, we repeated the modelling for these two epochs while restricting the secondary to a maximum brightness of 100% of the primary. We then repeated the astrometric fit, and found it to be incompatible with any physical orbit when considered with the other points. This means that the secondary’s increase in relative H band brightness to outshine the primary at the epochs 16 December 2019 and 28 December 2020 is the only solution that agrees with the Keplerian orbit of the objects.

Since these two points are at very similar orbital phases (0.88 and 0.91, respectively), we argue that this variability could be a result of the dynamical interaction of the secondary with the circumbinary disc, possibly through variable extinction, where the line-of-sight extinction changes towards one of the stars due to rearrangements in the circumbinary disc. In addition, our point-source flux estimate (f₂/f₁) might contain emission contributions from circumstellar gas or dust, in particular as the binary separation is still comparable to our interferometric beam size. Therefore, the observed f₂/f₁ brightness increase might correspond to an increase in excess emission near the location of the secondary as it approaches periastron, either through an accretion burst onto the secondary or from a hot spot at the inner edge of the circumbinary disc caused by the additional heating from the secondary. Indeed, smoothed particle hydrodynamics simulations of young, eccentric close binaries with circumbinary discs, such as those by Dunhill et al. (2015, for HD 104237) and Muzerolle et al. (2019, for DQ Tau), have suggested that dynamical interactions can cause a differential rate of accretion depending on the orbital phase, with the highest accretion peak occurring during the 10 − 20% of the orbit preceding periastron. The observed brightening, around orbital phase ∼0.9, could be consistent with this prediction.

6.3. Nature of the line-emitting region

As shown in Sect. 5, the geometry in the K band emission lines is substantially different from what we see in the continuum. The best fit to the circumstellar environment in the He I and Br γ lines indicates a strongly emitting Keplerian disc around the primary component. In this section we present three possible interpretations for the physical origin of the line emission: emission from the accretion region, an ionised gas accretion disc channelling circumstellar material to the star, or a decretion disc tracing material from the star being lost through stellar winds.

6.3.1. Accretion onto the primary

Brackett-γ emission is a common marker of magnetospheric accretion in YSOs and is especially useful for accretion rate determination, since the luminosity of the line (L_Brγ) has been shown to be directly related to the accretion luminosity L_acc over a wide mass range from brown dwarfs to Herbig Ae stars (e.g., Muzerolle et al. 1998; Calvet et al. 2004). This relation follows a power law that is strongly dependent on stellar mass (Fairlamb et al. 2017).

For Herbig Ae objects, Donehew & Brittain (2011) derived the following relation:

$$\begin{array}{c}\hfill log(L_{\mathrm{acc}}/L_{\odot })=(0.9\pm 0.2)log(L_{\mathrm{Br}\gamma }/L_{\odot })+(3.3\pm 0.7),\end{array}$$(5)

with similar values obtained by Mendigutía et al. (2011). From the above relation, the accretion rate Ṁ can be derived:

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\frac{G{M}_{\ast }\dot{M}}{R_{\ast }},\end{array}$$(6)

where M_* and R_* are respectively the mass and radius of the star. Using the Br γ luminosity for MWC 166 A calculated by Donehew & Brittain (2011) of L_Brγ = (11 ± 5)×10⁻³ L_⊙, we can estimate the accretion luminosity and hence the mass accretion rate. This results in L_acc ∼ 35 ± 65 L_⊙ and subsequently Ṁ ∼ (3.9±7.4) × 10⁻⁷ M_⊙ yr⁻¹, which is in line with the typical mass accretion rate of 2 × 10⁻⁷ M_⊙ yr⁻¹ found by Mendigutía et al. (2011).

However, it must be noted that, with a mass of 12 M_⊙, MWC 166 Aa is a Be star, leaving it outside of the regime where the L_Brγ – L_acc relation has been calibrated. Indeed, Eq. (5) provides a systematic overestimate of L_acc for Be stars, according to Donehew & Brittain (2011). Furthermore, the large uncertainties inherent to the estimate of L_Brγ by Donehew & Brittain (2011) result in very large errors on the propagated quantities. As a result, the above value for Ṁ should be taken to be an order-of-magnitude estimate.

The strong Br γ signal is consistent with accretion, as is the youth of the system. However, the geometry of the line emission does not appear to be consistent with direct stellar accretion. The radius of the emission originates from a region farther from the stellar surface than the ∼5R_⋆ that would be expected from accretion in Herbig Ae and Be stars (Bouvier et al. 2020). Furthermore, the geometry of the line emission is also stable over all epochs, with the red and blue lobes consistently located north-west and south-east of MWC 166 Aa, and with velocities consistent with Keplerian rotation, which is not expected of material being funnelled onto a stellar surface (Bouvier et al. 2007). As a result, the line emission does not appear to be tracing direct stellar accretion, but rather a process that is more spatially extended.

6.3.2. Inner gas accretion disc

The line emission might trace ionised gas in the inner region of the circumprimary gas disc. Based on hydrodynamic simulations, we expect that a stable circumprimary gas disc can exist out to one-third of the binary separation (Artymowicz & Lubow 1994). This disc could accommodate the mass transport from the large-scale mass reservoir seen in far-infrared excess emission, over a circumbinary disc to the star. For MWC 166 A, the upper limit on the radius is ∼0.9 au or ∼15 R₁, which is in agreement with the observed location of the Br γ emission (∼11.5 R₁). These discs are often found around Herbig Ae and Be systems, which would be consistent with the age of the system (see Sect. 6.4). The emission would then be tracing ionised gas in Keplerian rotation around the star, rather than direct accretion onto the star. Kraus et al. (2012) found an example of such a disc around another young B-type close binary system, V921 Sco, using similar techniques. The Br γ line profiles were found to be similarly strong and narrow to those presented in this study for MWC 166 Aa. They also were similarly variable in intensity and wing strength, while still being consistent with a Keplerian disc.

However, the lack of evidence of direct accretion in our observations means that this scenario cannot be confirmed. Further observations at shorter wavelengths may show evidence of accretion, for example by showing emission closer to the stellar surface.

6.3.3. Be decretion disc

As a third scenario, we consider that the observed circumprimary emission might be associated with a decretion disc, as observed around classical Be stars. Classical Be stars are defined as non-supergiant B stars that have shown Balmer line emission at some point in time (Collins et al. 1987). Generally very fast rotators, they typically host gaseous circumstellar discs characterised by a viscous decretion model (e.g., Lee et al. 1991; Carciofi 2011), with material ejected by radiation pressure from the star forming a Keplerian disc. This disc is sustained by periodic outbursts from the star (Grundstrom et al. 2011), which leads to variability on a range of timescales from days to years (Labadie-Bartz et al. 2017).

Variability in the emission line profile is also common in decretion discs, with violet-to-red (V/R) variations meaning the redshifted and blueshifted wings of the line are commonly found at different relative strengths at different epochs (Rivinius et al. 2013) – Fig. 3 shows that such variations are indeed present in the spectrum of MWC 166 A, and the results in Cols. (3) and (4) of Table 7 quantify this. In the simplest of Keplerian models, the flux intensities of each wing should be equal. The observed deviation from this profile is due to the flux intensity of the disc being azimuthally asymmetrical (Porter & Rivinius 2003), which is most likely due to temperature or gas density enhancements in specific regions of the disc. Our observations, comprising only four epochs, are insufficient to characterise the observed variations in V/R ratio as being periodic or quasi-periodic, but previous studies (Hanuschik et al. 1995) indicate that, in general, these variations are cyclic in the long term. Hummel & Hanuschik (1997) modelled a one-armed density wave precessing around the star as an explanation for these variations.

The rapidly changing V/R ratio is characteristic of Classical Be stars observed at intermediate inclinations (Catanzaro 2013). Such variations are generally not found in Herbig Be objects. Furthermore, the line-emitting region around the primary is rather extended, which is more consistent with decretion than accretion. We would assume Br γ emission from magnetospheric accretion to originate no more than 5R_⋆ from the star (Bouvier et al. 2020). However, we find that at all epochs, it is located well beyond the stellar surface. We estimated the radius of the emission by averaging the angular displacement in the blue and red wings of each line, with the results shown in Col. (5) of Table 7. These values do agree with each other in almost all circumstances, despite some slight variation between the two lines and epochs, with a general separation from the primary of ∼0.25 mas, or ∼11R₁. Comparing this value to other spectro-interferometric Be studies shows that it is consistent with that of a Classical Be disc, with the VLTI/AMBER survey by Cochetti et al. (2019) showing Br γ radii ranging between 3 − 13R_⋆. The lack of dust in Be discs also correlates with our findings from the continuum fit, where we found no evidence for circumprimary dust.

In light of the above points, the decretion disc model cannot be excluded for the circumprimary emission. Previous photometric studies of MWC 166 A found evidence for periodic dissipation and regeneration of the circumstellar envelope around the primary star over a period of months (Pogodin et al. 2006), which is consistent with models and observations of Be stars and adds weight to this interpretation of the line emission.

This would be in contrast to previous characterisations of MWC 166 Aa as a Herbig Be star (Fairlamb et al. 2015), and is also supported by our isochrone interpolation: from Fig. 9, it can be seen that the large mass of the primary means that it is evolved beyond the protostellar stage at all age estimates. This is further supported by the relatively small NIR excess in the spectral energy distribution of MWC 166 A, suggesting that the circumstellar environment has been at least partially cleared of material. The observed far-infrared excess of the system at ∼100 μm, however, is more supportive of the accretion disc scenario outlined in Sect. 6.3.2.

As a result, we can characterise the K band line emission as originating in a Keplerian or quasi-Keplerian circumprimary gas disc, as opposed to direct stellar accretion. Whether this is an accretion disc or a decretion disc is to be confirmed by future observations.

6.4. Constraining the age of the system

As Table 5 and Fig. 9 show, the primary star has already reached the main-sequence for all plausible age estimates, while the secondary is still in its pre-main-sequence stage at the youngest age estimates. By comparing the luminosities of the two stars at different ages to the flux ratios found from our interferometric observations, we can place some constraints on the age of the system.

We first performed the process described in Sect. 4.4 for every age isochrone in our CMD table between 100 kyr and 10 Myr. We then used the derived parameters to generate spectral energy distributions for both stars, using model atmospheres from Kurucz (1993). We then calculated values of f₂/f₁ at the central wavelengths of both the H and K bands, 1.65 μm and 2.2 μm, respectively, to enable a comparison with the corresponding values from our continuum interferometry (see Fig. 11).

Fig. 11. Flux ratio f2/f1 plotted against log(Age), for all isochrones (red lines), compared to interferometric flux ratios at the same wavelength (blue lines). The upper panel shows the flux ratio at 1.65 μm, while the lower panel shows the flux ratio at 2.2 μm – the central wavelengths for the H and K bands, respectively. The lighter-shaded areas illustrate the uncertainty on each quantity.

The H band flux ratio measurements show significant scatter between epochs, which might be due to scattered light contributions that are are important towards shorter wavelengths. As described in Sect. 6.2, PIONIER will likely be more strongly affected by these contributions, while MIRC-X and the long CHARA baselines should be able to resolve out extended flux and separate the components more reliably. Therefore, we adopt the MIRC-X value of f₂/f₁ = 0.484 ± 0.014 as our most reliable H-band flux ratio measurements, and mark this value with the blue line in the lower panel of Fig. 11.

The K band continuum GRAVITY data shows much less variability, and therefore we chose to take the weighted average of the GRAVITY flux ratios is plotted on the lower panel of Fig. 11, corresponding to an average flux ratio of f₂/f₁ = 0.399 ± 0.021.

Comparing the measured H and K band flux ratio with the model values, we can exclude ages for the system of ≲500 kyr. In the K band, this is well beyond the 3σ significance level. For the H band, this is only 2σ, but the range of ages < 500 kyr for which this is the case is very small, and only due to our relatively large mass uncertainties – which are themselves large because of the uncertainty on the system’s distance.

At such young ages, the NIR brightness of the secondary would exceed the primary, due to the low temperature of the secondary component, which is not supported by our interferometric observations. The upper age of the system is less well constrained, as the stars evolve very slowly once they reach the main-sequence (see Fig. 9). This can be seen in Fig. 11, where the models predict consistently f₂/f₁ ∼ 0.2 for ages beyond 1 Myr. Our average continuum K band flux ratio is larger than this, suggesting that the secondary component of the system is still in its pre-main-sequence stage. From our combined H and K band constraints, we estimate the system age to (7 ± 2)×10⁵ yr. This young age is also consistent with the presence of circumstellar material (see Sect. 6.1). Such an age is also consistent with either the accretion or decretion disc models described in Sects. 6.3.2 and 6.3.3.

7. Conclusions

In this study we have presented GRAVITY, PIONIER, and MIRC-X observations of MWC 166 A that resolve the system in the NIR H band and K band with milliarcsecond resolution. We derived the astrometry of the system at 13 epochs and calculated a first fully three-dimensional orbital solution for the system. This orbit differs substantially from the RV-only orbits of Corporon & Lagrange (1999) and Pogodin et al. (2006), having a period twice as long. We have subsequently constrained the dynamical system mass (17.1 ± 2.7 M_⊙ for the photometric distance of 990 ± 50 pc found by Kaltcheva & Hilditch 2000) and the distance of the system, with our results excluding previous parallax measurements of the distance where d < 500 pc. Furthermore, we have calculated, for the first time, the individual masses of the primary and secondary components of the system, which we find to be M₁ = (12.19 ± 2.18) M_⊙ and M₂ = (4.90 ± 0.52) M_⊙, respectively. We also find estimates for the other fundamental stellar parameters based on quadratic isochrone interpolation.

Furthermore, we see evidence for circumstellar emission, both in the dust continuum and in the He I and Br γ spectral lines, although they have different geometries. The geometry of the extended emission in the continuum is not well constrained, with the best fit corresponding to an over-resolved background halo. The variability in the continuum emission between epochs may be an indication of physical variability in the quantity of circumstellar dust, or it might indicate that the geometry is more complex than assumed in our model. Characterising it in more detail will require additional interferometric observations, ideally at mid-infrared wavelengths.

On the other hand, the geometry of the He I and Br γ line emission is well constrained by our observations. Our models show emission in these lines to be localised around the primary, where the redshifted and blueshifted wings are spatially displaced, consistent with gas in a circumprimary disc. The large spatial extent of the line-emitting regions (11.5R₁ for Brγ, 10.5R₁ for He I) and stable PA orientation are inconsistent with an origin in magnetospheric accretion or boundary-layer accretion, but support the hypothesis that the line emission is tracing an ionised gas disc. This gas disc might either be fed by mass infall from outside the binary or represent a decretion disc forming through mass loss from the primary.

Finally, we constrain the age of the system to (7 ± 2)×10⁵ yr, based on the measured flux ratio of the components. We find that the primary is a main-sequence Be star, while the secondary is a Herbig Be object still in the process of gravitational contraction onto the main sequence.

Acknowledgments

We acknowledge support from an STFC studentship (No. 109106G), an STFC Consolidated Grant (ST/V000721/1), and European Research Council Starting Grant “ImagePlanetFormDiscs” (Grant Agreement No. 639889) and Consolidated Grant “Gaia-BIFROST” (Grant Agreement No. 101003096). Travel support was provided by STFC PATT grant ST/S005293/1. MIRC-X received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 639889). JDM acknowledges funding for the development of MIRC-X (NASA-XRP NNX16AD43G, NSF-AST 1909165) and MYSTIC (NSF-ATI 1506540, NSF-AST 1909165). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in Ochsenbein et al. (2000). This research has made use of the Jean-Marie Mariotti Center JSDC catalogue available at http://www.jmmc.fr/catalogue_jsdc.htm. This research has made use of the Jean-Marie Mariotti Center OiDB service available at http://oidb.jmmc.fr. Analysis in this work was based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme(s) 098.C-0910(A) (GRAVITY); 190.C-0963(A), 190.C-0963(B), 102.C-0701(B), 104.C-0737(A), 104.C-0737(B), 104.C-0737(C), 106.21JU.001, 106.21JU.002, 106.21JU.003 (PIONIER). This research has made use of the PIONIER data reduction package of the Jean-Marie Mariotti Center available at http://www.jmmc.fr/pionier. 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. We thank the referee for their insightful comments.

References

Appendix A: PMOIRED line models

Below we show the full set of He I and Br γ line models for the four GRAVITY epochs, discussed in Sects. 3.2 and 5. The figures comprise nine panels, laid out equivalently to those in Fig. 10.

Fig. A.1. 14 March 2017, He I

Fig. A.2. 14 March 2017, Br γ.

Fig. A.3. 27 April 2017, He I

Fig. A.4. 27 April 2017, Br γ

Fig. A.5. 11 January 2018, He I

Fig. A.6. 11 January 2018, Br γ

Fig. A.7. 6 February 2018, He I

Fig. A.8. 6 February 2018, Br γ

All Tables

All Figures

	Fig. 1. Visibilities (and associated residuals) of MIRC-X models. At V ≳ 0.6, the purely point-source model overshoots the observed data points substantially.
In the text

	Fig. 2. Closure phases (and associated residuals) of MIRC-X models.
In the text

	Fig. 3. Continuum-normalised spectra around the He I and Brγ lines. The different epochs have been offset for clarity. Epoch-dependent variations are visible. Dates are given in the format YYYY-MM-DD.
In the text

	Fig. 4. Observed continuum visibilities (black) and corresponding models (red) plotted against spatial frequency for all epochs. The model fit includes extended background emission. Dates are given in the format YYYY-MM-DD.
In the text

	Fig. 5. Observed continuum closure phases (black) and corresponding models (red) plotted against spatial frequency for all epochs. The model fit includes extended background emission. The scaling on the vertical axis was adjusted for each epoch. Dates are given in the format YYYY-MM-DD.
In the text

	Fig. 6. Corner plot showing the possible correlations between free parameters (ρ, θ, F2 = 100 ⋅ f2/ftot and F3 = 100 ⋅ fext/ftot, respectively) for epoch 6 February 2018 of the continuum GRAVITY data, where the model fit includes extended background emission.
In the text

	Fig. 7. Astrometric orbit solutions derived using both the ORBITX code (blue line) and the grid-search code (red line). The primary star is kept fixed at the origin, and the x and y axes show displacement in right ascension and declination, respectively. The dotted lines connect the ascending and descending nodes of each orbit.
In the text

	Fig. 8. Radial velocity measurements of MWC 166 A taken in the period 1994–2005 plotted against orbital phase. The blue and red fitted curves correspond to the orbits specified in Cols. (3) and (4) of Table 4, respectively. The dotted black line shows the velocity of the system’s centre of mass (V0) for the grid-search method, and the solid grey line shows V0 for the ORBITX orbit.
In the text

	Fig. 9. Hertzsprung-Russell diagram showing PARSEC 1.2S isochrone tracks for each age. Superimposed are the interpolated L, Teff values for the primary star (large circles) and secondary star (triangles), for all ages. The black cross shows the location of the primary’s parameters as determined by Fairlamb et al. (2015).
In the text

Fig. 10. Results of the Br γ line modelling, for the epoch 14 March 2017. (a): (u, v)-coverage for the observations associated with the epoch, coloured by baseline pair. (b): telluric-corrected normalised flux (labelled NFLUX, black lines), overplotted with flux in the best-fit model (red line). (c)–(e): data from each GRAVITY exposure (black lines), overplotted with quantities computed from best-fit model (red lines). The observables are closure phase for each telescope triplet (T3PHI), differential phase for each baseline (DPHI), and visibility for each baseline (|V|), respectively. (f)–(h): brightness distribution corresponding to the best fit, for three representative wavelengths. (i): synthetic line strengths and profiles for the two spectral components in the model (red and blue), in addition to the monochromatic continuum flux associated with the primary component (silver) and the secondary component (gold). The parameters corresponding to the fits can be found in Table 6.

In the text

	Fig. 11. Flux ratio f2/f1 plotted against log(Age), for all isochrones (red lines), compared to interferometric flux ratios at the same wavelength (blue lines). The upper panel shows the flux ratio at 1.65 μm, while the lower panel shows the flux ratio at 2.2 μm – the central wavelengths for the H and K bands, respectively. The lighter-shaded areas illustrate the uncertainty on each quantity.
In the text

	Fig. A.1. 14 March 2017, He I
In the text

	Fig. A.2. 14 March 2017, Br γ.
In the text

	Fig. A.3. 27 April 2017, He I
In the text

	Fig. A.4. 27 April 2017, Br γ
In the text

	Fig. A.5. 11 January 2018, He I
In the text

	Fig. A.6. 11 January 2018, Br γ
In the text

	Fig. A.7. 6 February 2018, He I
In the text

	Fig. A.8. 6 February 2018, Br γ
In the text