{'fraction_non_alphanumeric': 0.030675072138800566, 'fraction_numerical': 0.020309140031980002, 'mean_word_length': 4.436650679456434, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ARMADA I: Triple Companions Detected in B-Type Binaries alpha Del and nu Gem

Ground-based optical long-baseline interferometry has the power to measure the orbits of close binary systems at ~10 micro-arcsecond precision. This precision makes it possible to detect"wobbles"in the binary motion due to the gravitational pull from additional short period companions. We started the ARrangement for Micro-Arcsecond Differential Astrometry (ARMADA) survey with the MIRC-X instrument at the CHARA array for the purpose of detecting giant planets and stellar companions orbiting individual stars in binary systems. We describe our observations for the survey, and introduce the wavelength calibration scheme that delivers precision at the tens of micro-arcseconds level for<0.2 arcsecond binaries. We test our instrument performance on a known triple system kappa Peg, and show that our survey is delivering a factor of 10 better precision than previous similar surveys. We present astrometric detections of tertiary components to two B-type binaries: a 30-day companion to alpha Del, and a 50-day companion to nu Gem. We also collected radial velocity data for alpha Del with the Tennessee State University Automated Spectroscopic Telescope at Fairborn Observatory. We are able to measure the orbits and masses of all three components in these systems. We find that the previously published RV orbit for the inner pair of nu Gem is not consistent with our visual orbit. The precision achieved for these orbits suggests that our ARMADA survey will be successful at discovering new compact triple systems to A/B-type binary systems, leading to better statistics of hierarchical system architectures and formation history.


INTRODUCTION
Binary systems are of vital importance for understanding the physical properties of stars by providing the opportunity to directly measure stellar masses. Precision visual binary orbits can be combined with single or double-lined radial velocity (RV) data to solve for the physical parameters of the system (i.e., masses, physical orbital elements, mutual inclinations for tertiaries, position on HR diagram). Along with the importance of providing mass and stellar evolution information, visual binary orbits at the micro-arcsecond (µas) level of precision can open up the possibility to search for additional inner companions to the outer binary system. When a binary orbit is monitored over a long enough time scale, unseen inner companions will impart an additional "wobble" motion on top of the binary motion as the inner pair orbits its center-of-mass. Once instrumental precision reaches the ∼10µas level, this method can be used to detect companions down to the regime of Jupiter-mass planets. The power of this method opens up exciting new areas of parameter space in which to search for exoplanets and stellar companions.
For mass ratios >0.1 for the inner pair, tertiary occurrence rate is expected to be roughly 30% for A/B-type binaries  Table 2 for details of the calibrators used. STS is our internal calibration source, implemented in 2019May.  (Bonneau et al. 2006) The MIRC-X combiner measures visibilities, differential phase, and closure phase of our targets. Normally, one employs frequent observations of nearby calibrator stars to measure visibility loss due to time-variable factors such as atmospheric coherence time, vibrations, differential dispersion, and birefringence in the beam train. For ARMADA our main interest is differential astrometry between two components of a binary system within the interferometric field-of-view (both components generally unresolved). Since closure phase is immune to atmospheric effects, and extra dispersion in differential phase can be fit with a polynomial, we are able to observe for ARMADA without the use of traditional calibrators by fitting to closure and differential phase. This is crucial for ARMADA operations, as we are able to spend most of the night observing targets rather than spending time on calibrator sources. We still use sparse on-sky calibrators for our wavelength calibration (described in section 2.3) as well as for a rough calibration of visibilities. In May 2019, the MIRC-X team implemented an internal light source called Six Telescope Simulator (STS, described in Anugu et al. 2020). This higher SNR source is now used for wavelength calibration instead of the on-sky calibrator in dates following its implementation. We used the MIRC-X data pipeline (version 1.3.3) to produce OIFITS files for each night, described in Anugu et al. (2020). This pipeline and its documentation is maintained on Gitlab 1 . These nights were reduced with the "spectral-differential" method of the MIRC-X pipeline for computing differential phase. This method first removes the group delay from the raw phase, and then computes differential phase as the phase(i+1) -phase(i) where i, i+1 are neighboring wavelength channels. We reduced most of our data with the number of coherent integration frames (ncoh) of 10, oifits max integration time of 60 seconds, and bispectrum bias correction applied. Since the companion to κ Peg is near the edge of the interferometric field-of-view (>150 mas for latter epochs), its binary phase signal is varying faster on some baselines than the 60-second integration time. Hence we reduce the κ Peg epochs from 2019 June onward with an oifits max integration time of 10 seconds. Since this object is bright, the high SNR allows us to combine fewer frames into a measurement of phase. Data taken in 2018 September showed signs of a vibration present within our combiner, with a quick loss of coherence reported by the pipeline. This forced us to reduce any 2018 September data with a lower number of ncoh (3 frames, instead of 10).

Fitting Binary Star Differential Astrometry
For each MIRC-X night we fit to the following binary model of complex visibility, V : The free parameters for this binary model include a uniform disk for the primary and secondary to form visibilities V 1 and V 2 ; a binary separation in right ascension (R.A.) and declination (DEC) -(α, δ); a monochromatic flux ratio between the two components f ; as well as a bandwidth smearing parameter b = 1/R, where R is the resolution of the disperser and Γ = sinc [b(uα + vδ)]. The location on the uv-plane is denoted by parameters u and v.
Since we do not use the standard CAL-SCI sequence of observing, our squared visibilities are poorly calibrated. We thus use the closure phase and differential phase observables to fit our binary positions for each epoch. To inform our uniform diameter (UD) values in our fits, we choose a well calibrated epoch for each target (i.e., a dataset that is taken near to the nightly calibrator star). This choice has little affect on the differential astrometry, since these targets are all mostly unresolved. For α Del we use the 2019Jul29 epoch (calibrator HD 176437) to fit for uniform diameters. We see no χ 2 improvement in the fits by letting the diameters vary, and the two companions of this triple system prefer diameter values close to point sources. Hence we fix all three diameters to U D = 0.5 mas, which is near the resolution limit in H-band. Gordon et al. (2019) measured a uniform diameter of 0.407 ± 0.022 mas for the A component of α Del, consistent with our fixed values. ν Gem also showed no improvement in the fits by letting the UDs vary for the night of 2017Sep30 (calibrator HD219080). Hence we again fix all three diameter to U D = 0.5 mas. For the well-calibrated 2019Aug06 epoch of κ Peg (calibrator HD 886), we measure diameters of U D 1 = 0.730 ± 0.002 mas and U D 2 = 0.733 ± 0.003 mas. We then fix the uniform diameters for all κ Peg epochs to U D = 0.7 mas. The flux ratio values we report are the fitted values from these well-calibrated epochs.
In the case of α Del and some epochs of ν Gem we also detect flux from a third component in the system, as judged by significant residuals in a fit of the binary model to the observables. In this case, an extra component is added to our complex visibility model: The symbols have their same meaning as in the binary model, except now in our notation f 1 , f 2 , f 3 are the flux contributions from each component and f 1 + f 2 + f 3 = 1. Note there are now two differential positions in R.A. and DEC -the primary to the secondary (α 12 , δ 12 ), and the primary to the tertiary (α 13 , δ 13 ). To find our best differential astrometry solution on a given night, we first perform a wide grid search in R.A. and declination with step sizes of 0.1 milli-arcseconds to find the minimum χ 2 solution. We then perform a non-linear least squares fit using the lmf it package in Python to narrow in on the best solution (Newville et al. 2016). To search for a third component, we perform an additional grid search of the tertiary component and fit for the expected binary (initial guesses informed by our first grid search) on each point of this grid ( Figure 2). After obtaining the approximate location of all three components, we perform a global fit with our triple complex visibility model starting from the best guesses from the grid searches. Figure 1 shows an example fit for one of our epochs of α Del.
HD196867 Closure Phase HD196867 Differential Phase Figure 1. We show an example fit of our triple model to our observables closure and differential phase of α Del from 2019Jul30. The two panels show our best fit (dashed red line with crosses) to closure phase (left) and differential phase (right). The dashed grey line shows the best-fit binary model, which is a significantly worse fit to our data than the triple model. Each square in the plots represents one of MIRC-X's 20 closing triangles or 15 baselines (the 6 telescopes are designated as E1-E2-W1-W2-S1-S2). This is a single 60-second measurement, of which we typically have ∼10 per observation of a target. To search for flux from triple companions, we first perform a grid search for the brighter binary (left). Due to the signal from a tertiary companion, this initial χ 2 map is quite messy. We then do a grid search for a third component, with a fit for the binary companion performed at each iteration (right). The location of the binary component is denoted by the star in the plot. Including this third component vastly improves the chi-square maps. Once we narrow in on the best solutions for each component on the grid, we perform a joint least-squares fit to refine the solutions.

Binary Grid Search Tertiary Grid Search
We convert our astrometry solutions from differential R.A. and DEC to a separation and position angle east of north (ρ, θ). After finding the best-fit ρ and θ, we estimate errors by performing a 2D scan in R.A. and declination around the best solution. The error ellipses we report are then a fit to the region bound by χ 2 = χ 2 best ± 1. Since these error ellipse sizes are conservative estimates, we further scale these errors down to raise χ 2 red = 1 after performing our full set of orbital fits described in section 3.

Wavelength Calibration
In most previous astrometric experiments (e.g., PTI-PHASES: Colavita et al. 1999;Muterspaugh et al. 2010, Keck-ASTRA: Woillez et al. 2010, VLTI-PRIMA: Sahlmann et al. 2013, VLTI-GRAVITY: Gravity Collaboration et al. 2017), a differential delay line with laser metrology is used to measure the pathlength difference between broadband fringe packets from the two stars under study. For ARMADA, we use a spectrometer with spectral resolution of R = λ ∆λ ∼ 200 so that each star's interferogram is about 200 fringes long. This allows the packets from the stars to overlap each other, creating distinctive wavelength dependent variations in the fringe phase. These variations can be directly related to the separation of the two stars and do not require any kind of differential delay line. The downside is that we must know the effective wavelength of each spectral channel to high precision. MIRC-X is limited to a wavelength precision knowledge of ∆λ λ ∼ 10 −3 (Monnier et al. 2012). For a 100 mas binary, this would impose a limitation up to 100 µarcsec astrometry, too large for our ultimate goals of detecting the astrometric signal from orbiting exoplanets.
We employ an extra calibration step to bring our wavelength precision knowledge to the 10 −4 level. Our team has utilized a custom-built 6-beam optical etalon system to calibrate our astrometry each ARMADA night. Each etalon consists of a thin, ∼2mm-thick piece of glass with parallel sides and 50% reflective coatings on both sides. By having an etalon in each beam being a slightly different thickness (∆x = 6 ± 2µm), our combiner detects multiple interference packets that look similar to a binary star. We measure the spectrometer wavelength using the same data pipeline and methodology as we are using for the science targets. We designed and built a special, thermally-stable etalon holder which is placed into the 6 CHARA beams each night. Our model for the etalon signal is similar to that of a binary signal in Equation 1, though the binary separation is now replaced by the differential thicknesses of the etalons: In this equation n λ is the index of refraction for infrasil 301: n λ = [0.6961663λ 2 /(λ 2 −0.0684043 2 )+0.4079426λ 2 /(λ 2 − 0.1162414 2 ) + 0.8974794λ 2 /(λ 2 − 9.896161 2 ) + 1] 1/2 , ∆s is the differential thickness between two etalons of different beams, and f j is the fraction of flux in a given reflection. V 1 is the primary beam visibility, while now the bandwidth smearing parameter is defined as Γ = sinc[b(n λ ∆s/λ)]. The etalons create an infinite number of reflections which can be modeled in theory, though we find that including more than two terms (j=1 and j=2) is not necessary for improving our χ 2 fits. Figure 3 shows an example dataset for a measurement of etalon signal from the internal STS source.
A systematic error in the wavelength solution for a given night will affect the measured value of binary separation. We thus use our etalon calibration data to bring each night to the same astrometric scale. Each ARMADA night, we take etalon data to generate an astrometric correction factor to apply to the measured binary separations. Before the implementation of the internal STS source, these data were taken on-sky with the calibrator star. We fit for differential thickness between the 6 etalons for each beam (free parameters are five ∆s values between one reference beam, 2nd order polynomial for differential phase dispersion on each baseline, and the monochromatic flux ratio lost in a reflection). We then choose a high SNR reference night for ARMADA, and scale each night's binary separations based off this reference. The choice of the ARMADA reference night is arbitrary, since we are dealing only with internal astrometric consistency between ARMADA nights. In future work, we plan to use a shared binary source between MIRC-X and VLTI-GRAVITY in order to carry out an absolute calibration of MIRC-X wavelengths (GRAVITY achieves precise absolute wavelength calibration with an internal fourier transform spectrometer source). Currently, there is an estimated 0.25% wavelength precision for MIRC-X (Monnier et al. 2012). This systematic applies to the measured binary separations and needs to be taken into account when combining our presented astrometry with high precision astrometric data from other instruments.
We calculate a single "etalon correction factor" for each night based off of the slope of the 15 REFERENCE vs NIGHT etalon optical path differences (OPDs -one per baseline). Since our etalons are in a thermally-stable holder, any change in our etalon model comes from the MIRC-X wavelength solution. This scaling is then applied to the separation ρ of the binaries for the night, and hence is done as the final step after finding the best differential astrometry solution described in the previous section. The separations that we report have our etalon calibration applied. In Fig 3 we show how this scale factor changes with time across our ARMADA datasets. The maximum effect that this correction has on our separations is a 4e-3 factor (400 micro-arcseconds for a 100mas binary), with the median scaling being at a 7.8e-4 factor (78 micro-arcseconds for a 100mas binary). The nights which have a correction of ∼4e-3 are worse than the expected systematic for MIRC from Monnier et al. (2012). The likely culprit of this high correction factor is a change in the detector readout windowing, which was set in 2018Feb to avoid a bad pixel in the fringe window. Since the spectrograph has a slightly different optical magnification across the field of view, this change of position of the fringe window in the detector would affect the fringe spatial frequency and thus the wavelength solution. This "new" fringe window was unchanged for ARMADA data taken between 2018Feb -2018Aug, which is consistent with the higher etalon correction values. The error bars on the etalon correction factor are computed by bootstrapping the etalon datasets in time, and taking the standard deviations of the resulting scale factors computed. The error bars on our etalon factors are all smaller than the reported errors for the separation from astrometry, and hence do not add error to the astrometric solution.
To test that our etalon wavelength correction scheme is working, we fit binary orbits for ARMADA both with and without the correction applied. Since our data only record a fraction of the outer orbital period, we also include historical data from the Washington Double Star (WDS) Catalog (Mason et al. 2001). Table 3 shows how the mean and median residuals to these orbit fits change (the full orbit fitting routine is presented in section 3). The median residual to the best fit orbit decreases significantly for all three targets after applying our etalon calibration. Though the median residual of ν Gem decreases significantly, it shows the least improvement in mean residual when applying the etalon correction (due to a few higher residual points). This added noise is possibly due to time-varying resolved structure in the Be-disk of this system, as discussed in Section 5.2. Figure 4 shows the residual fit in R.A. and DEC for κ Peg, before and after the etalon correction is applied.  We show the R.A. and DEC residuals of κ Peg both before and after applying our etalon correction factor. The mean and median residual values both decrease, bringing us to our desired ∼10 micro-arcsecond precision which is needed to eventually detect substellar companions. The most significant correction in this case is a single outlier point which was >100 micro-arcseconds off from the expected orbit before the etalon factor is applied.

ORBIT FITTING
Once we have our measured binary separations and position angles for each night, we are able to fit a Keplerian orbit to the data. Since we are aiming for high precision differential astrometry, we need to account for the precession of North when combining position angles measured by MIRC-X to historical data in the WDS catalog. The MIRC-X pipeline already defines the uv-plane in ICRS using the Python "astropy.coordinates" package (Astropy Collaboration et al. 2013, thus accounting for Earth orientation. Note that we also make extensive use of the Numpy package in our Python routines (Harris et al. 2020). We correct the position angles of the WDS data to a common J2000 reference with P.
where t is the year of observation (Siregar & Kuncarayakti 2010). The Campbell elements (ω, Ω, e, i, a, T , P ) describe the Keplerian motion of one star of a binary system relative to the other. Those symbols have their usual meanings where ω is the longitude of the periastron, Ω is the position angle of the ascending node, e is the eccentricity, i is the orbital inclination, a is angular separation, T is a time of periastron passage, and P is orbital period. For near-circular orbits T and ω become ill-defined, adding seemingly large errors to these orbital parameters as the two parameters are correlated. In these cases we also report a time of maximum radial velocity T rv,max , which is more tightly constrained. When including RV data, we also fit to the semi-amplitudes K and system velocity γ. The longitude of periastron ω is traditionally reported for the secondary when fitting to visual binary orbits alone. The convention when combining RV orbits is to report ω of the primary, which is flipped by 180 • . We have RV data for both κ Peg and α Del, and hence we report the ω of the primary for these orbits. For ν Gem we only present the visual orbit, and so we report ω of the secondary (noted in the table of orbital elements). For visual orbits, there is a 180 degree ambiguity between ω and Ω. Our RV information breaks this degeneracy for κ Peg and α Del, and we report the Ω < 180 • for ν Gem.
For nonlinear least-squares fitting, we use the Thiele-Innes elements to describe our Keplerian orbits. As described in Wright & Howard (2009), these elements convert (ω, Ω, i, a) to linear parameters (A, B, F, G). When fitting a system of three components, we assume the 3-body system is hierarchical with the wide companion orbiting the center-of-mass of the inner pair. This means that our orbit model is simply a sum of the outer + inner Keplerian orbits. Since the outer orbits we present in this paper are significantly larger than the inner orbits (>200 times larger in orbital period), this hierarchical model is a reasonable assumption. When flux from the third component is not detected, our orbital elements are then describing the "wobble" motion of one star about the center-of-mass of the inner orbit. In this case the angular semi-major axis a wob of the tertiary component describes the size of the wobble motion, where one would need to know the mass ratio to figure out the true angular semi-major axis of the inner pair. In the case where we detect flux from all three components, we then modify our orbital elements to also include the mass ratio since we are then able to measure a 1 and a 2 of the inner semi-major axis (a inner = a 1 + a 2 ). We then do a joint fit to the outer orbit, the wobble motion, and the inner visual orbit.
We again use the Python lmf it package for non-linear least-squares fitting of our data (Newville et al. 2016). To constrain the outer binary orbits, we include historical data from WDS. Based on the large scatter of these data about their best fits, we assign circular errors of radius 10 milli-arcseconds to the WDS data. We use the ORB6 catalog for initial guesses of our orbital parameters for the outer pair. Once we find the best fit for the outer binary, we begin searching for the inner companion. To do so, we vary the inner orbital period and fit circular orbits to each fixed period as shown in Figure 5. Once the best inner period is detected, we refine our search further by performing a joint outer + inner fit with all orbital parameters varying.
α Del ν Gem Error bars for the fitted orbital parameters are normally estimated in lmf it from the covariance matrix, but since the orbital elements P and e are nonlinear we instead determine posterior distributions on our orbital parameters with a Markov chain Monte Carlo (MCMC) fitting routine. We carry out MCMC fitting using the Python package emcee developed by Foreman-Mackey et al. (2013). We use our best-fit orbital elements as a starting point for our 2*N params walkers, where the starting point for each walker is perturbed about its best fit value. We assume uniform priors on all of our orbital elements. The quoted error bars on our orbital elements are the standard deviations of the posterior distributions, and corner plots of the posteriors for the inner and outer orbits show correlations between parameters.
4. KAPPA PEG: VERIFYING ARMADA ASTROMETRY Before presenting our new astrometric detections in systems α Del and ν Gem, we wanted to test our ARMADA calibration scheme on a well-known system to verify our precision astrometry. Muterspaugh et al. (2006) published a high precision orbit on the triple star system κ Peg (HD 206901, HIP 107354, WDSJ21446+2539) with their interferometric PHASES survey. This system consists of a wide ∼0.2" pair of F5 subgiants, and an inner component around the brighter star with a 6-day orbital period. We observed this tertiary system over the course of a year with our ARMADA survey, and compare our inner orbit to that obtained by Muterspaugh et al. (2006). For historical reasons, the brighter component in H-band was designated the B component of the outer A-B system -a nomenclature that the previous work kept. To stay consistent with most previous work on κ Peg and make a direct comparison of orbital elements, we will keep this designation when reporting orbital elements. Hence the brighter star is the Ba+Bb inner subsystem.
We obtained 8 new data points for this system throughout 2018-2019, taken after the MIRC-X optics upgrade of 2018Sep. In Figure 6, we show the MIRC-X and PHASES data on our best fit to the outer AB orbit of ∼11.5 years. We first ignore the PHASES data for a comparison, fitting only to the MIRC-X data and historical WDS data -important for constraining the outer orbit. Our high precision astrometry captures the motion of the brighter component about the center-of-mass of its inner pair. We show our fit to this "wobble" motion after subtracting out the binary motion, in order to visually compare the astrometry between the PHASES and MIRC-X datasets. We present our best-fit MIRC-X orbital elements in Table 5, along with the previously published elements for comparison. Figure 7 is a good demonstration of the improvement in interferometric data quality over the past decade between PHASES and ARMADA. Our median residual to the binary+wobble motion is an amazing 6.32 micro-arcseconds after we apply our etalon wavelength calibration, when fitting to ARMADA+WDS datasets. This demonstrated precision is about a factor of 10 improvement over PHASES, which is very promising for our ultimate goal with the ARMADA survey to detect circumstellar giant planets in A/B-type binary systems.
We also combine our new astrometric datasets with the PHASES data and radial velocity data presented in Muterspaugh et al. (2006) for a full combined fit. Our update to the best orbit of κ Peg is presented in the final column of Table 5. Note that Muterspaugh et al. (2008) reported updated orbital elements for this system, as there was a sign flip present in the original analysis. We use these updated orbital elements for the PHASES+RV column of this table. Since we view the "wobble" motion of the bright component, along with the single-line RV motion of this component, we are able to compute orbital parallax for the system. Our value agrees well with the Hipparcos distance of 34.2 ± 0.9 pc (van Leeuwen 2007), as well as the previously published value for PHASES. We also update the masses of all three components in the system. Figure 9 shows the corner plots of our posterior distributions from the MCMC routine. We see no significant correlations between the parameters of the inner orbit with the parameters of the outer orbit -hence we split these corner plots into outer and inner orbital elements to increase clarity. In general, the outer long-period orbit shows greater correlations between the orbital parameters. This is not surprising, as the coverage is not as complete for the ∼11.5 year period as compared to the inner 6-day period. For near circular orbits, the correlation between the angle of periastron and the time of passage through periastron is expected -as these parameters become ill-defined. This is especially apparent in the inner orbit. There is also a strong correlation between the inclination and semi-major axis for the outer orbit, which is normal for partial coverage of a visual orbit.
We were not able to detect the flux from the Bb component of the system. Our measured mass for this component of 0.814 ± 0.046M could be consistent with that of a white dwarf remnant of a massive star, though Muterspaugh et al. (2006) claim to see evidence of a third set of lines in their Keck-HIRES spectra (although RVs of this third set were not measured in that work). This makes it more likely that the Bb component is a late G or early K dwarf, which should imply a flux ratio of ∼7% with respect to the Ba component in the near infrared. MIRC-X is easily capable of detecting companions of this flux ratio, though the signal from the near 1:1 wide binary dominates in this case. Our rough visibility calibration scheme for the ARMADA survey also increases the difficulty of detecting faint nearby companions (expected semi-major axis is ∼2.5 mas). Though we found some tentative detections of the third companion in a few of our nights, they all appeared to be related to residual structure in the χ 2 maps. These tentative detections also did not align with the expected position angle from the "wobble" orbit, hence we concluded that they were not real. Figure 6. κ Peg is a triple star followed by both ARMADA and PHASES surveys. We show the data from MIRC-X and PHASES together, along with the best fit to the outer orbit. We also plot historical data obtained from the WDS catalog. The zoomed inset to a portion of our ARMADA data gives an idea of the size of our error ellipses. . We show our best-fit orbit to the inner "wobble" of κ Peg due to the tertiary companion, after subtracting out the outer binary motion. We compare the results of the PHASES best-fit orbit from Muterspaugh et al. (2006) (left) to our MIRC-X+WDS fit (right). Over the ∼year time baseline of following this object with MIRC-X we are achieving a median residual level of 6.32 micro-arcseconds, which is very promising for our ultimate goal of detecting giant planets in binary systems with the ARMADA survey. Here the RV motion of the bright primary due to the inner companion is subtracted out. We also show data of the brighter component phased to the 5.9 day inner period (lower right). Figure 9. We show corner plots of the MCMC posterior distributions produced from emcee for the outer (left) and inner (right) orbital elements for κ Peg. Since the outer long-period orbit is less well-constrained than the short-period inner orbit, there are more obvious correlations between the parameters. Since the inner orbit is near circular, there is a strong correlation between ω2 and T2.  Muterspaugh et al. (2008). We compute these from the elements given in that paper. c Flux ratios in K are derived from Keck AO imaging in Muterspaugh et al. (2006), rather than PHASES data.

α Del
The star α Del (HD196867, HIP101958, WDSJ20396+1555) is bright (V = 3.8, H = 3.9) with a spectral type of B9IV (Morgan & Keenan 1973). α Del was first determined to have a binary companion by Wickes (1975). That study used the Mt. Wilson 60-inch telescope to find a ∼200 milli-arcsecond binary companion with a magnitude difference of about 2 in the visual passband. McAlister (1977) then followed up on this new system with speckle interferometry on the Mayall 4-meter telescope, confirming the new detection and in subsequent studies adding additional points to the binary orbit. The binary orbit has since been well constrained with speckle data, with WDS reporting nearly 100 data points from 1974-2015. The Sixth Catalog of Orbits of Visual Binary Stars (ORB6, Malkov et al. 2012) lists a 17 year period for the binary orbit, with an eccentricity of 0.47. The ORB6 orbit is graded a 2 on a scale of 1-5, meaning most of the revolution is covered and the orbit quality is considered good.
We report the discovery of a third object in this system. We find that the B component is itself a short period subsystem, with an orbit of 30 days. This additional "wobble" can clearly be seen with our astrometric precision of the outer binary orbit in Figure 10. We also detect the flux from the newly discovered Bb component in our MIRC-X data, which allows us to see the motion of both Ba and Bb components in the new pair about their center-of-mass. With 11 epochs, we are able to fit a full triple orbit model. To help constrain the outer orbit, we also include 70 epochs of historical data from WDS.
To better constrain the new inner orbit and obtain a measurement of orbital parallax, we collected new RV data from the Tennessee State University 2 m Automated Spectroscopic Telescope (AST) and its echelle spectrograph at the Fairborn Observatory (Eaton & Williamson 2007). See Gardner et al. (2018) for a discussion of the spectroscopic observations from this telescope and their velocity reduction. From 2020 January through October we obtained 55 RVs, covering the orbit of the inner pair. We identified 2 sets of lines from the spectra -a broad-line component (v sin i ∼145 km s −1 ) for the primary detected with a line list for A stars, and a second very narrow component (v sin i ∼7 km s −1 ), which can be seen with both an A star line list and a solar star line list. The first RV component is stationary, and thus we attribute it to the primary A-component. The second RV component shows the ∼30-day motion of the Ba-Bb pair, which can be used to compute a 1 sin i of the inner orbit. It is not immediately obvious whether the measured RV is that of the Ba or Bb component. Since it is detected in both the A-star and solar line lists, this suggests a component which has a late A or early F spectral class. Our orbital parallax and mass determinations in Table 10 are consistent with this picture if the measured RVs are that of the less massive Bb component. The measured masses in this scenario would imply a mid-to-early A star for the Ba component. A rapidly rotating A-star would make this component difficult to detect in the composite spectra, and would explain why this third component is not easily visible despite it being the more massive star of the inner pair. A close inspection of the mean composite line profile of the stars, which in total is only about 1 percent deep and was obtained with the A-star line list, shows a very weak, broad asymmetry, with a v sin i of roughly 60 km s −1 that shifts in the opposite velocity direction to the component that we have identified as Bb. This provides supporting evidence for our conclusion that the narrow lines are from component Bb.
Since we know the inclination, the angular semi-major axis of Ba+Bb, and the motion of Ba about the center-of-mass from astrometry, and the a 2 sin i of Bb from RV, we are able to measure a 2 of Bb in both physical and angular units. This gives us a measurement of orbital parallax which can be used to compute masses of all three components without depending on an outside measurement of distance to the system. Our distance of 78.5 ± 2.2 pc is consistent with the Hipparcos value of 77.8 ± 2.7 pc (van Leeuwen 2007). The Gaia DR2 parallax of 66.8 ± 2.4 pc (Gaia Collaboration 2018) does not agree well with our measured orbit, though this is likely because DR2 does not yet take multiplicity into account when computing the parallax. We measure three masses of 3.83 ± 0.33 M , 1.82 ± 0.15 M , and 1.49 ± 0.12 M for the primary, Ba component, and Bb component respectively. We report the orbital elements for the outer pair and this newly detected tertiary companion in Table 9. Table 10 gives the physical properties of these stars which can be deduced from the single-line RV + visual orbits. Figure 12 shows posterior distributions of the outer and inner orbital elements. As for κ Peg in the previous sections, there are no significant correlations between the parameters of the inner orbit with the parameters of the outer orbit. We again split the parameters up in the plot for clarity. The inner orbit of α Del is particularly well constrained. Since the inner orbit it near circular, there is the expected correlation between T and ω.    The motion of Ba+Bb about their center-of-mass allows us to work out the mass ratio of this inner pair, which in turn leads to measurements of all three masses in the system. Figure 12. We show corner plots for the inner (left) and outer (right) orbital elements of α Del. The inner orbit is particularly well constrained in this case, with little correlation among orbital parameters.

ν Gem
ν Gem (HD 45542, HIP 30883, WDSJ06290+2013) is a bright (V =4.14, H=4.43) binary with a B-type primary and a Be star for the visual secondary at ∼0.1". This B+Be pair has a visual magnitude difference of about 1, with a period of about 19 years reported in the ORB6 catalog. Jarad et al. (1989) discovered a third star in this system with a period of 40-days, though the authors pointed out that they are not confident in this period value. Rivinius et al. (2006) confirmed this inner spectroscopic binary from radial velocity variations seen between 1994-2003, with a best-fit period of 53.72 days. These authors note that the circumstellar Balmer emission and shell lines do not take part in the 53.75-day period, concluding that the shell star is component B. The hierarchy of this system adopted from Rivinius et al. (2006) is Aa+Ab for the inner pair, with the B component of the wide AB orbit being the Be star (note that this notation differs from the WDS catalog, which denotes the wide pair presented here as Aa-Ab with a B component at >100 arcsec). We report the first astrometric detection of the inner Aa+Ab orbit as a "wobble" to the outer binary motion, with a period in agreement with Rivinius et al. (2006).
Recovering the flux from the inner component is difficult for this system since the signal is dominated by the bright wide binary, and our group delay tracker was noisier than usual for these datasets. When coupled with poor calibration, this makes a detection of a close (<1 mas in some epochs) companion difficult. We were able to confidently detect the inner companion in two of our MIRC-X epochs near apastron. To add additional data points to the inner visual orbit, we used 2 old MIRC datasets from 2015 and 2016. These datasets were taken in the lower resolution PRISM mode (R∼42), meaning that the wide component is outside the interferometric FOV on most baselines. The stable fringe tracking on these nights and traditional TARGET-CALIBRATOR calibration sequence make it easier to find the inner component. The MIRC datasets were reduced with the standard MIRC data pipeline in IDL described in previous MIRC papers (e.g. Monnier et al. 2012 Schaefer et al. 2010). On 2016Nov14 we used calibrators 71 Ori, HD23183, and HD37329. We scaled the wavelengths produced by the pipeline by a factor of 1.004 to bring to an absolute scale as suggested by Monnier et al. (2012). For these old MIRC datasets, we used an IDL binary grid search routine 2 modified to fit the system as a triple (Schaefer et al. 2016). The component diameters were fixed as point sources. Although the wide binary was outside the interferometric field-of-view on many baselines, we still fit for this wide component in our MIRC data. However, since the old MIRC epochs do not have our etalon wavelength calibration data, we are unable to bring these data to the same scale as our more recent MIRC-X nights. We assign the MIRC epochs errors of 0.4 mas (consistent with the largest etalon scale factor correction of 0.4% = 0.4 mas for a 100 mas binary). At separations of a few milli-arcseconds for the inner pair, any wavelength scalings at the 1e-3 level have no effect on the measured separation within reported error bars. The MIRC positions are reported along with the MIRC-X astrometry in Tables 11 and 12. In order to measure orbital parallax, we initially combined our astrometric data with the single-line RV provided by Rivinius et al. (2006). Our joint visual+RV orbit yielded a "wobble" semi-major axis of 1.47 mas, inclination of 81 degrees for the inner orbit, inner period of near 54 days, and a RV semi-amplitude of 35 km s −1 for the primary. These values lead to a distance of 125 pc, implying masses of the system much too low to be consistent with B-type stars. With an outer semi-major and period of 83 mas and 6974 days, the total mass in the system is only ∼3 M at a distance of 125 pc. In private communication with the authors we learned that there are new RV data for this system, and that an upcoming paper will update the RV orbit. The published RVs in (Rivinius et al. 2006) were measured assuming an SB1 orbit, when in fact the spectral lines are a blend of the two inner components, leading to inaccurate RVs. Klement et al. (in prep) will include a full analysis of the 3 stars in this system, combining our astrometric new data with their updated RV analysis of all three components. Hence we present the visual orbit of ν Gem, without joint fitting of the problematic RVs. In Table 14 we assume the Hipparcos distance of 167 ± 8 pc (van Leeuwen 2007) when computing masses.
When fitting to the WDS + ARMADA data for ν Gem, we find that the best-fit orbital parameters change between fitting the inner "wobble" alone vs. fitting the inner visual orbit. As can be seen in Table 13, our "wobble" orbit prefers a more eccentric solution. Our median residual to the outer binary + inner wobble motion for this orbit is 14.8 µarcsec. When coupling this motion with the position of the Ab component, however, our median residual increases to 40 µarcsec and the solution becomes less eccentric. This added residual in the wobble could imply motion from an additional companion in the system (which can be fit out by driving up eccentricity). In fact, there is accumulating evidence that nearly all Be stars are close binaries (e.g. Klement et al. 2017Klement et al. , 2019. The Hipparcos distance of this system leads to a mass sum which is already too low for B-type stars, implying that additional companions would need to be low mass -although new RV data might update this distance. Recent work suggests that the companions to Be stars are generally low-mass, stripped down cores of former mass donors (Wang et al. 2017;Bodensteiner et al. 2020). Though additional companions are still possible, more epochs are needed to confirm the orbit of a fourth body in this system. It is also possible that the residuals can be explained by resolved time-varying structures in the Be disk. Since the "wobble" depends on the measurement from star Aa -B, this structure would add residual motion to our measured "wobble" orbit.
Our measured masses for the three components of this system are 2.7 ± 0.4 M , 2.5 ± 0.4 M , and 1.8 ± 0.4 M for components Aa, Ab, and B respectively. The higher error bars come from the high error on parallax from Hipparcos. These masses assume that the Hipparcos distance for the system is accurate, which often is not the case for close binaries and triples. Rivinius et al. (2006) listed B6III and B8IIIe as the preferred spectral types for the primary and wide component, though the spectral classification of Be stars in general is very uncertain due to weak photospheric lines. Classifications for this system in literature have ranged from mid-to-late B-type stars. Our reported masses are too low for these spectral classifications, likely implying that the distance is too low. Gaia DR2 lists an even lower distance of 155 ± 14 pc (Gaia Collaboration 2018), though the lack of a multiple star solution makes this measurement untrustworthy. Klement et al. (in prep) will be able to measure an orbital parallax by combining their updated RV analysis with the astrometry presented in this work, which will lead to a more accurate measurement of masses.
The orbital elements of the outer pair of ν Gem have higher errors bars when compared to α Del. This is due to worse coverage of this ∼19 year period. Figure 15 shows posterior distributions of the outer and inner orbital elements. Again, there are no significant correlations between the parameters of the inner orbit with the parameters of the outer orbit -so we split them up for clarity. Looking at the outer elements, one can see many correlations between parameters due to incomplete coverage of this outer orbit. Better coverage of this wide orbit may alter the outer orbital elements, particularly constraining semi-major axis and inclination more confidently. These elements are important for constraining the total mass in the system, and better coverage of this wide orbit would lead to better mass constraints.
Trv,max (MJD) -58510.7 ± 0.2 58511.60 ± 0.34 58511.34 ± 0.09 a Since we do not include RV, there is a 180 • degeneracy. We report the Ω < 180 • . b ±0.25% from absolute wavelength precision c "wobble" semi-major Figure 13. We detected the inner tertiary component to the system of B-type binary ν Gem. We show the best fit to the ARMADA+WDS data for the outer binary pair -we zoom in to the ARMADA portion of the orbit to show that there is clearly an additional "wobble" to this binary motion. Here we are fitting only to the outer binary motion + wobble. (Center) We also detect flux from the Ab component in four MIRC-X epochs, and in two epochs of MIRC data. We plot the best fit orbit when fitting only to this Aa-Ab pair. Visually, one can see that the best-fit orbit does not agree well with the orbit fit of the wobble motion alone. The wobble motion prefers a higher eccentricity orbit than that derived from the visual data alone, with a significantly larger period. (Right) We show the combined fit, when coupling the wobble motion to the visual orbit. The median residual of the orbit increases, due to the residuals from the wobble motion. This extra residual is potentially due to time-varying resolved structure from the Be star disk, or additional companions in the system. Figure 15. We show corner plots for the inner (left) and outer (right) orbital elements of ν Gem. Since the outer 19-year period is not very well constrained, many of the orbital parameters show correlations with each other. Better coverage of the outer orbit will lead to tighter constraints. The inner orbit parameters show fewer correlations, with the strongest being between the angle and time of passage through periastron -expected for a near circular orbit. ainner (au) 0.48 ± 0.02 a The Hipparcos distance is assumed (van Leeuwen 2007). This affects the measurement of masses and physical semi-major axes. Klement et al. (in prep) will use our data to obtain a better measurement of orbital parallax.

Implications for formation
These triple systems exemplify the power provided by complete orbital information, particulary in the case of the well-constrained α Del orbit. Our analysis shows that the Ba-Bb binary orbit of α Del is retrograde with respect to its orbit about the primary, with a mutual inclination of 159.8 ± 0.5 • . Though not unique (see e.g. Tokovinin & Latham (2020), retrograde orbits are atypical. Unlike many known high inclination triples α Del is not subject to Kozai-Lidov oscillations (Lidov 1962;Kozai 1962). With near equal mass ratio and an eccentric outer orbit, the system stretches the limits of current secular theory (Naoz 2016). Thus we verify via n-body simulations using Rebound (Rein & Liu 2012) that the full orbital solution is stable against secular instability, and merely undergoes nodal precession and small amplitude eccentricity oscillations. Notably the current oscillating orbital elements are very close to the median values found over 1000s of orbits.
That the system is retrograde, stable against large inclination excitation, and compact, suggests a possibly violent dynamical interaction in the past. Formation within the same disk via fragmentation is strongly disfavored given the mutual inclination and short orbital periods. (Kratter et al. 2010;Hall et al. 2017). Formation via turbulent fragmentation on large scales (or a combination of both modes of fragmentation) requires substantial migration either via gas-drag or dynamical interactions; the former case would favor smaller mutual inclinations (Lee et al. 2019). In the latter case, the tightness of the orbits favors the ejection of one or more previously bound objects to serve as the requisite sink of energy and angular momentum.
An expected consequence of dynamical instability is the misalignment of orbital and stellar spin vectors. While the obliquities of the three components of α Del cannot be measured directly, we can infer from measured rotational velocities that alignment is not favored for two of the three components. For Ba, the measured v sin i ≈ 60km/s is consistent with typical spin rates for mid-A stars if it is aligned with either the Ba-Bb orbit (v ≈ 160 km/s at 22 • ) or the BaBb-A orbit (v ≈ 180 km/s at 161 • ) (Zorec & Royer 2012). However, with v sin i ≈ 7 km/s, the Bb component would have to be a very slow rotator (v ≈ 20 km/s) if it shared Ba's obliquity, unusual for a star above the Kraft break. Similarly, the primary component's v sin i ≈ 145km/s is hard, though not impossible, to reconcile with an inclination of 161 • : a nearly pole on orbit with v ≈ 445 km/s implies that the star is rotating near breakup velocity. With our current data, it is not yet possible to rule out the expected overluminosity or distortion consistent with such rapid rotation.
Without the inclusion of RV information, ν Gem is harder to constrain from a formation history perspective due to the ambiguity of the position angles of the ascending node (Ω) for the inner and outer orbits. In contrast with α Del, our ν Gem data is consistent with nearly co-planar orbits -though RV data will break the degeneracy. Co-planar orbits often have a formation history with sequential epochs of disk fragmentation followed by migration (Tokovinin & Moe 2020). Disk instability becomes increasingly likely at higher stellar masses (Kratter et al. 2008). The upcoming RV analysis for ν Gem, including measurements of v sin i for some components, will increase the clarity of mutual inclinations and potential stellar obliquities in this system (Klement et al., in prep).

SUMMARY AND FUTURE WORK
We started the ARMADA survey with the MIRC-X instrument at the CHARA array with the goal of detecting giant planets on ∼au orbits and previously unseen low-mass companions orbiting individual stars of binary systems. In this paper, we introduced the observational methods and calibration scheme which can bring us to ∼10 micro-arcsecond precision for our differential binary orbits. Our newly implemented etalon module for precise wavelength calibration was demonstrated to improve systematic errors in binary separation on the triple star test system κ Peg. We achieved a median precision level of 6.3 micro-arcseconds for data taken after our 2018Sep optics upgrade of MIRC-X, and demonstrated a vast improvement in astrometric performance compared to previous similar interferometric surveys. We presented the updated RV+visual orbit of this triple star system. Though our astrometric precision is about a factor of 10 better than previous work, there is still potential to do better in approaching the fundamental limits of interferometric observations. Ireland & Woillez (2018) predict submicroarcsecond atmospheric and shot-noise limited precision for the ∼100 milli-arcsecond binaries presented in this paper. Hence another factor of ∼10 improvement is theoretically possible, but would require calibration of pupil registration at the <1 mm level projected onto primary telescope mirror space. Future ARMADA papers will included more detailed studies in the factors currently limiting precision.
Confident detections of giant planets will take a longer time-baseline and higher number of epochs, but are beginning to see the "wobble" signature from the gravitational tug of previously unseen short-period tertiary companions. Since these companions are stellar, we often detect the flux from the stars as well and measure the mass ratio of the inner pair. Combined with single-line RV data, we are able to measure orbital parallax along with the masses of all three components in these systems. We made the first discovery and measured the orbital elements of a 30-day companion to the B-type binary α Del, which includes our new MIRC-X/ARMADA data as well as new RV points from the Fairborn Observatory. We also detected for the first time the inner visual orbit of B-type triple star ν Gem. This system is of particular importance, since it will become part of an ongoing study to understand multiplicity in Be star systems. Our detection of the inner orbit is crucial for solving the physical parameters of the system, and an upcoming paper by Klement et al. (in prep) will characterize this system by including new RV data to this inner and outer orbit.
Our data is consistent with α Del being born of a violent dynamical instability / ejection event, while being more uncertain for ν Gem without RV data. The rich dynamics of these systems demonstrate the power of full orbit solutions for revealing the origins of triple systems (both planetary and stellar). Future discoveries that contain, for example, more compact inner binaries may also prove useful for constraining basic tidal evolution models.
As is evident from these first three systems, our ARMADA survey will be extremely efficient at detecting <1 au stellar companions in wide binary systems. In future work, we plan to publish our full list of newly detected triple systems with the ARMADA survey. When combined with RV of the outer and inner pairs, we will be able to fully characterize these orbits, including masses (which can be measured using the methods of this paper), along with mutual inclinations (if combined with RV), and evolution history on an HR diagram.
The residuals to our best-fit orbits are also promising for detecting ∼au brown dwarfs and giant planets in interesting regimes. We target mainly A/B-type stars for ARMADA, where the ∼au giant planet occurrence is difficult to measure due to weak and broad spectral lines of these stars. Binary systems themselves are crucial regimes for searching for planets, as circumstellar planets in close binaries are difficult to detect. Finding planets in binary systems provides clues to the formation timescales and channels for giant planets. With more epochs, our ARMADA survey should be able to detect giant planets in these regimes. Facilities: CHARA, Fairborn Observatory