# ξ Tauri: a unique laboratory to study the dynamic interaction in a compact hierarchical quadruple system

J. Nemravová, P. Harmanec, M. Brož, D. Vokrouhlický, D. Mourard, C.

Hummel, C. Cameron, J. Matthews, C. Bolton, H. Božić, et al.

► To cite this version:

J. Nemravová, P. Harmanec, M. Brož, D. Vokrouhlický, D. Mourard, et al. "ξ Tauri: a unique laboratory to study the dynamic interaction in a compact hierarchical quadruple system. Astronomy and Astrophysics - A&A, 2016, 594, pp.A55. 10.1051/0004-6361/201628860". hal-03473851

HAL Id: hal-03473851

https://hal.science/hal-03473851

Submitted on 13 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ξ Tauri: a unique laboratory to study the dynamic interaction in a compact hierarchical quadruple system☆,☆☆

J. A. Nemravová¹,*,*,*, P. Harmanec¹, M. Brož¹, D. Vokrouhlický¹, D. Mourard², C. A. Hummel³, C. Cameron⁴,
J. M. Matthews⁵, C. T. Bolton⁶, H. Božić⁷, R. Chini⁸,⁹, T. Dembsky⁸, S. Engle¹⁰, C. Farrington¹¹, J. H. Grunhut³,
D. B. Guenther¹², E. F. Guinan¹⁰, D. Korčáková¹, P. Koubský¹³, R. Kříček¹, R. Kuschnig¹⁴, P. Mayer¹,
G. P. McCook¹⁰, A. F. J. Moffat¹⁵, N. Nardetto², A. Prša¹⁰, J. Ribeiro¹⁶, J. Rowe¹⁷, S. Rucinski⁶, P. Škoda¹³,
M. Šlechta¹³, I. Tallon-Bosc¹⁸, V. Votruba¹³, W. W. Weiss¹⁴, M. Wolf¹, P. Zasche¹, and R. T. Zavala¹⁹

*(Affiliations can be found after the references)*

Received 5 May 2016 / Accepted 27 May 2016

## ABSTRACT

**Context.** Compact hierarchical systems are important because the effects caused by the dynamical interaction among its members occur on a human timescale. These interactions play a role in the formation of close binaries through Kozai cycles with tides. One such system is ξ Tauri: it has three hierarchical orbits: 7.14 d (eclipsing components Aa, Ab), 145 d (components Aa+Ab, B), and 51 yr (components Aa+Ab+B, C).

**Aims.** We aim to obtain physical properties of the system and to study the dynamical interaction between its components.

**Methods.** Our analysis is based on a large series of spectroscopic photometric (including space-borne) observations and long-baseline optical and infrared spectro-interferometric observations. We used two approaches to infer the system properties: a set of observation-specific models, where all components have elliptical trajectories, and an *N*-body model, which computes the trajectory of each component by integrating Newton's equations of motion.

**Results.** The triple subsystem exhibits clear signs of dynamical interaction. The most pronounced are the advance of the apsidal line and eclipse-timing variations. We determined the geometry of all three orbits using both observation-specific and *N*-body models. The latter correctly accounted for observed effects of the dynamical interaction, predicted cyclic variations of orbital inclinations, and determined the sense of motion of all orbits. Using perturbation theory, we demonstrate that prominent secular and periodic dynamical effects are explainable with a quadrupole interaction. We constrained the basic properties of all components, especially of members of the inner triple subsystem and detected rapid low-amplitude light variations that we attribute to co-rotating surface structures of component B. We also estimated the radius of component B. Properties of component C remain uncertain because of its low relative luminosity. We provide an independent estimate of the distance to the system.

**Conclusions.** The accuracy and consistency of our results make ξ Tau an excellent test bed for models of formation and evolution of hierarchical systems.

**Key words.** binaries: close – binaries: spectroscopic – binaries: eclipsing – stars: kinematics and dynamics – stars: fundamental parameters – supernovae: individual: ξ Tauri

## 1. Introduction

Binaries and multiple systems play a crucial role in our understanding of the formation, stability, and evolution of stars and their hierarchies, starting from simple binaries up to galaxies.

Of all known binaries, those that eclipse have represented the most useful group because until recently, an accurate determination of component masses and radii was possible primarily for them. For binaries with components of different masses, a common origin of the system also provided a stringent test of the models of stellar evolution. At the same time, however, this fact represented an unpleasant selection effect, especially for binaries

with hot components and rapid rotation: we have only observed them roughly equator-on so far.

The recent rapid advances in optical interferometry allowing the usage of longer baselines, co-phasing of more telescopes, and longer integration times provide the opportunity of obtaining accurate basic physical properties for non-eclipsing binaries as well. It is possible to obtain the spatial orbit of these binaries and derive their accurate orbital inclination. In combination with radial-velocity (RV) curves, this allows determining component masses and the absolute value of the semi-major axis. Since the interferometric orbit provides the angular value of the semi-major axis, we also obtain an estimate of the distance of the binary that is completely independent of the photometric distance modulus. In the most favourable cases, long-baseline interferometry can also provide independent estimates of the component radii.

Many binaries are members of multiple systems (Eggleton & Tokovinin 2008). When it is possible to derive masses of more than two components, not only the nuclear but also the dynamical evolution of such systems can be studied. It has been suggested that the formation of triple systems, containing a compact binary accompanied by a distant component,

<footnote>* Full Tables D.1–D.7 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/594/A55</footnote>

<footnote>** Based on data from the MOST satellite, a former Canadian Space Agency mission, jointly operated by Microsatellite Systems Canada Inc. (MSCI; formerly Dynacon Inc.), the University of Toronto Institute for Aerospace Studies and the University of British Columbia, with the assistance of the University of Vienna.</footnote>

<footnote>*** Coresponding author: J. A. Nemravová, e-mail: jana.nemravova@gmail.com</footnote>
was dynamically very exciting. During the evolution, gravitational interactions of the three stars are expected to excite the eccentricity of the binary through the Kozai mechanism, which brings them close to each other. Later, tides stabilise the system by preventing the Kozai-pumped eccentricity from further increasing and revert the trend to circularisation (e.g. Eggleton & Kiseleva-Eggleton 2001; Fabrycky & Tremaine 2007). Even though we cannot observe the systems at their dynamically violent youth, we can still appreciate some degree of dynamical evolution produced by continued gravitational interactions of the three stars. To compare predictions of the theory with observations, the mutual orientation of orbits with respect to each other is required, that is, their inclinations and the longitudes of ascending nodes. These are available only for objects for which an astrometric orbit is known. This in turn can only be obtained with interferometry.

We here investigate one such system, the unique and rare close quadruple system $\xi$ Tau, whose favourable orbital geometry as well as the luminosity ratios between its components allow determining physical properties of the system and its components with high precision. Possible dynamical effects in the system can be studied as well. $\xi$ Tau (2 Tau, HD 21364, HIP 16083, and HR 1038) is a hierarchical quadruple system consisting of two sharp-lined A stars that undergo binary eclipses, a more distant broad-lined B star, and a much more distant F star. The visual magnitude $V = 3.72$ mag, the declination of $9^{\circ}44'$, and the quite accurate HIPPARCOS parallax $15.6 \pm 1.04$ mas (van Leeuwen 2007) make $\xi$ Tau an easy and interesting target for a wide range of instruments and observational techniques.

The binary nature of the system was discovered by Campbell (1909). The wide orbit was first resolved by Mason et al. (1999) through speckle interferometry. All later available speckle-interferometric observations were analysed by Rica Romero (2010), who derived an astrometric orbit. The inner triple system was first mentioned by Fekel (1981), who quoted orbital periods of $7.15$ d and $145.0$ d based on a private communication from C. T. Bolton. The orbital elements of the triple subsystem were published in a catalogue by Tokovinin (1997). More accurate elements were given in a preliminary report by Bolton & Grunhut (2007), who obtained periods of $7.1466440(49)$ d and $145.1317(40)$ d. They were also the first to note that the inner binary is an eclipsing system, based on HIPPARCOS photometry. Hummel et al. (2013) reported a solution of the $145.2$ d orbit based on interferometric observations. The first detailed, but still preliminary study of $\xi$ Tau was published by Nemravová et al. (2013). These authors analysed numerous spectral, photometric and interferometric observations and discovered the apsidal motion of the $145.2$ d orbit with a period $224 \pm 147$ yr. They were able to separate the spectra of the two A stars and the broad-lined B star.

The system is quite complex, hence we briefly summarise its orbital elements and the properties of its components based on our analysis as presented in following sections in Table 1. It serves only to introduce the system and is not to be confused with our results.

This paper represents a comprehensive study of the system, based on analyses of a huge and unique body of spectral, photometric, and spectro-interferometric and astrometric data. Each type of observation is first analysed separately by standard means (Sects. 3–6), and the results are then critically compared in Sect. 7. Using them as the initial starting point, we then present the N-body model of the whole quadruple system, in which the mutual interactions of the orbits are also modelled. This is a new approach that tries to embrace almost all available

<table><caption>**Table 1.** Brief summary of orbital elements and properties of components of $\xi$ Tau.</caption><thead><tr><th rowspan="2">Quantity</th><th rowspan="2"></th><th colspan="3">Orbit</th></tr><tr><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>P</td><td>(d)</td><td>7.14664</td><td>145.12</td><td>18 630</td></tr><tr><td>e</td><td></td><td>≤0.01</td><td>0.21</td><td>0.564</td></tr><tr><td>i</td><td>(deg)</td><td>86.8</td><td>86.6</td><td>-24.4</td></tr><tr><td>a</td><td>(R<sub>☉</sub>)</td><td>25.3</td><td>233</td><td>6097</td></tr><tr><td>Cpts.</td><td></td><td>Aa+Ab</td><td>A<sup>1</sup>+B</td><td>AB<sup>2</sup>+C</td></tr><tr><td>Sp.T.</td><td></td><td>B9V+B9V</td><td>+B5V</td><td>+F5V<sup>3</sup></td></tr><tr><td>m</td><td>(M<sub>☉</sub>)</td><td>2.25+2.13</td><td>+3.73</td><td>+0.92<sup>4</sup></td></tr><tr><td>V</td><td>(mag)</td><td>5.46+5.53</td><td>+4.25</td><td>+7.55<sup>3</sup></td></tr></tbody></table>

**Notes.** This table serves only for introductory purposes and does not present our final results. P denotes the orbital period, e the eccentricity, i the inclination, Cpts. are the components of an orbit, Sp.T. the spectral type, m the mass, and V the apparent magnitude in the Johnson V filter. (1) A denotes the inner eclipsing system Aa+Ab. (2) AB denotes the intermediate system A+B. (3) Based on the magnitude difference listed in the HIPPARCOS and Tycho catalogue (ESA 1997). (4) The evident mismatch between the spectral type and its mass is due to the high uncertainty of both parameters.

pieces of information and provides the best description of the geometry and dynamics of the system to date (see Sect. 8). Finally, we recall some results of a simple perturbation theory in Sect. 9, which allows us to understand the principal dynamical effects revealed by the numerical model in Sect. 8.

We denote the individual components and orbits of the system as follows: components Aa and Ab are the primary and secondary of the close eclipsing subsystem revolving in a $7.15$ d orbit, labelled 1. Component B is the broad-lined star of spectral type B, revolving with the close pair in the $145$ d orbit, labelled 2. Finally, we denote the faint and very distant F-type star as component C and its $51$ yr orbit with the triple subsystem as orbit 3.

## 2. Observations and reductions

Here we provide only basic information about the observational material at our disposal. More details on the datasets and their reductions are provided in Appendices A–C.

Throughout this paper we use a shortened form of heliocentric Julian dates, reduced Julian dates given as $RJD = HJD - 2400 000.0$.

### 2.1. Spectral observations

The series of spectroscopic observations that has previously been used by Nemravová et al. (2013) was complemented with more recent ones secured at Ondřejov and La Silla: they were made with the echelle spectrograph FEROS (Kaufer et al. 1999), and at Cerro Armazones with the BESO spectrograph (Steiner et al. 2008; Fuhrmann et al. 2011). Four archival ELODIE echelle spectra were also used (Moultaka et al. 2004). With this rich collection of electronic spectra, we no longer needed the early RVs from the David Dunlap Observatory (DDO) photographic spectra that were used by Nemravová et al. (2013). The spectra were primarily used to obtain RV measurements of all three components of the close triple subsystem. The journal of all available spectra with the number of measured RVs for the components of the inner triple subsystem is listed in Table 2. More details on the spectra and their reductions can be found in Appendix A.
<table><caption>**Table 2.** Journal of spectroscopic observations.</caption><thead><tr><th>ΔT<br>(RJD)</th><th>N<br>Aa/Ab/B</th><th>Δλ<br>(Å)</th><th>R</th><th>Ins.</th></tr></thead><tbody><tr><td>49 300.7–52 670.5</td><td>37/37/37</td><td>4357–4568</td><td>10 800</td><td>DDO</td></tr><tr><td rowspan="3">51 960.3–53 637.6</td><td>04/04/04</td><td>4270–4523</td><td rowspan="3">42 000</td><td rowspan="3">ELO</td></tr><tr><td>04/04/04</td><td>4759–4991</td></tr><tr><td>04/04/04</td><td>6260–6735</td></tr><tr><td rowspan="3">55 041.9–55 867.6</td><td>13/13/13</td><td>4270–4523</td><td rowspan="3">48 000</td><td rowspan="3">BES</td></tr><tr><td>13/13/13</td><td>4759–4991</td></tr><tr><td>13/13/13</td><td>6260–6735</td></tr><tr><td>55 579.4–56 357.3</td><td>34/34/34</td><td>4270–4523</td><td>19 200</td><td>OND</td></tr><tr><td>56 579.4–56 889.6</td><td>05/04/05</td><td>4274–4508</td><td>19 200</td><td>OND</td></tr><tr><td>55 579.3–55 645.3</td><td>02/02/02</td><td>4378–4632</td><td>17 700</td><td>OND</td></tr><tr><td>55 579.3–56 357.3</td><td>20/20/20</td><td>4753–5005</td><td>19 300</td><td>OND</td></tr><tr><td>56 527.6–56 592.5</td><td>05/05/05</td><td>4759–4991</td><td>21 500</td><td>OND</td></tr><tr><td>56 527.6–56 889.6</td><td>14/14/14</td><td>6260–6735</td><td>14 000</td><td>OND</td></tr><tr><td>55 561.3–56 357.3</td><td>58/58/59</td><td>6255–6767</td><td>12 700</td><td>OND</td></tr><tr><td>55 597.4–55 980.3</td><td>19/19/22</td><td>6497–6688</td><td>14 000</td><td>LIS</td></tr><tr><td rowspan="3">56 555.7–56 564.7</td><td>12/12/12</td><td>4270–4523</td><td rowspan="3">48 000</td><td rowspan="3">FER</td></tr><tr><td>12/12/12</td><td>4759–4991</td></tr><tr><td>12/12/12</td><td>6260–6735</td></tr></tbody></table>

**Notes.** For each instrument, ΔT refers to the time span between the first and the last measurement, N gives the number of RVs measured for components Aa, Ab, and B, Δλ is the wavelength interval covered by the spectra in question, and R is the spectral resolution. In column “Ins.”: DDO – David Dunlap Observatory 1.9 m reflector, Cassegrain CCD spectrograph; ELO – Haute-Provence Observatory 1.2 m reflector, echelle ELODIE CCD spectrograph; BES – Cerro Amazones Hexapod Telescope, BESO echelle CCD spectrograph; OND – Ondřejov Observatory 2 m reflector, coudé CCD spectrograph; LIS – Lisbon Observatory of the Instituto Geográfico do Exército, reflector, CCD spectrograph; FER – La Silla 2.2 m reflector, Feros echelle CCD spectrograph.

Radial velocities measured on the available spectra (see Sect. 3.2) are listed in Table D.1.

## 2.2. Photometric observations

The photometry that has previously been used by Nemravová et al. (2013) was complemented by very accurate observations acquired almost continuously over two weeks with the MOST satellite (Walker et al. 2003) and by another series of Johnson UVB observations from Hvar. Additionally, we analysed the photometric minima published by Zasche et al. (2014).

The MOST satellite monitored ξ Tau over 16 days almost continuously. It acquired 21 525 observations that after the initial reduction by the MOST team were still affected by two systematic effects: the stray light from the Earth atmosphere, which introduced narrow peaks with separation ≈101 min; this is the MOST orbital period. The other effect was the relaxation time after the change of the observed field, during which the CCD had to reach thermal equilibrium. This manifests itself by a slowly decreasing offset that typically lasts several tens of minutes. The first effect was, with the exception of few observations during eclipses, removed with a low-passband Butterworth filter (Butterworth 1930). The second effect forced us to neglect all observations secured before RJD = 56 522. The remaining 18 510 observations were then analysed.

A journal of available photometric observations is listed in Table 3, and more details on the observations and data reductions can be found in Appendix B.

The reduced UVB photometric observations acquired at the Hvar Observatory, at the South African Astronomical

Observatory, the Four College APT, and photometric observations acquired with the MOST satellite are listed in Tables D.2–D.5.

## 2.3. Interferometric observations

The system was observed by four different spectro-interferometers: the Mark III Stellar Interferometer<fnref n="1" /> (Mark III) (Shao et al. 1988), the Navy Precision Optical Interferometer (NPOI) (Armstrong et al. 1998), the Visible spEctroGraph and polArimeter (VEGA) (Mourard et al. 2009) mounted at the Centre for High Angular Resolution Astronomy (CHARA) (ten Brummelaar et al. 2005), and the Astronomical Multi-BEam combineR (AMBER) (Petrov et al. 2007) attached to the Very Large Telescope Interferometer (VLTI) (Glindemann et al. 2004). A journal of the spectro-interferometric observations is listed in Table 4. The phase coverage of orbits 1 and 2 with all spectro-interferometric observations is shown in Fig. 1. Details on the spectro-interferometric observations and their reduction are provided in Appendix C.

Reduced spectro-interferometric observations from all four instruments are listed in the form of calibrated squared visibility moduli in Table D.6 and closure phases are provided in Table D.7.

# 3. Spectroscopy

The spectral lines of all three components of the triple subsystem (i.e. orbits 1 and 2) of ξ Tau are clearly seen in all available spectra. Component C was not detected in any of the spectra at our disposal because its relative luminosity is lower than 1%, which is beyond the detection limit of the available spectra. Attempts to detect its spectral lines were carried out through spectral disentangling and a comparison of the near-infrared spectra with synthetic profiles, both with null results.

Two different approaches to derive the orbital elements of the triple subsystem of ξ Tau were used. The first was a direct analysis of RVs measured with the method described in Sect. 3.2, and the second was the spectral disentangling (Simon & Sturm 1994; Hadrava 1995) in Sect. 3.3.

Additionally, we derived the basic radiative properties of ξ Tau using the comparison of the synthetic to observed and separated spectra (i.e. obtained through the spectral disentangling).

## 3.1. RVs measured by comparing the observed and synthetic line profiles

The RVs were derived using an automatic method based on the comparison of synthetic and observed spectra that searches for the best match with the optimisation of $\chi^2$ given by

$$ \chi^2 = \sum_{i=1}^{N_I} \left( \frac{I_{\text{OBS}}(\lambda_i) - \sum_{j=1}^{N_C} I_{\text{SYN},j}(\lambda_i, RV_j)}{\sigma_i} \right)^2, \quad (1) $$

where $I_{\text{OBS}}$ is the observed spectrum, $I_{\text{SYN},j}$ the synthetic spectrum of the jth component, $N_I$ is the number of discrete elements of the digitised spectrum, $N_C$ is the number of the components of the system, $RV_j$ is the radial velocity of the jth component, and $\sigma_i$ the standard deviation of the ith point of the observed spectrum, which was estimated from the continuum and adopted for the whole spectrum.

<footnote n="1">¹ Decommissioned in 1992.</footnote>
<table><caption>**Table 3.** Journal of photometric observations.</caption><thead><tr><th>N</th><th>ΔT<br>(RJD)</th><th>Passbands</th><th>Comp/Check</th><th>Instrument</th></tr></thead><tbody><tr><td>441/451/452</td><td>46 324.6–56 882.6<sup>1</sup></td><td>UBV</td><td>4 Tau/6 Tau</td><td>HVAR</td></tr><tr><td>69</td><td>47 909.6–48 695.0</td><td>V(H<sub>p</sub>)<sup>2</sup></td><td>all-sky</td><td>HIPP</td></tr><tr><td>26/26/26</td><td>55 569.3–55 579.4</td><td>UBV</td><td>6 Tau/4 Tau</td><td>SAAO</td></tr><tr><td>131/133/135</td><td>55 883.9–55 956.8</td><td>UBV</td><td>4 Tau/6 Tau</td><td>VILL</td></tr><tr><td>18 510</td><td>56 222.0–56 238.0</td><td>MO</td><td>all-sky</td><td>MOST</td></tr></tbody></table>

**Notes.** For each row, *N* is the number of observations in each of the filters used, Δ*T* is the time span covered by each dataset, column “Passbands” shows the photometric filters used, column “Comp/Check” lists the names of comparison and check star used. *UBV* denote the Johnson filters, and *MO* denotes the broad-band filter of the MOST satellite. (1) Only three observations were taken before RJD = 54 116, all at RJD = 46 324.

(2) The original HIPPARCOS *H*<sub>p</sub> broad-band observations were transformed into the Johnson V filter following Harmanec (1998). However, for the light-curve solutions the limb-darkening coefficients corresponding to the original HIPPARCOS passband were used. Instruments: HVAR – Hvar Observatory 0.65 m Cassegrain reflector, photoelectric photometer; HIPP – The ESA HIPPARCOS Astrometric Mission; SAAO – South African Astronomical Observatory 0.5 m Cassegrain reflector, Lucy photoelectric photometer; VILL – the Four College 0.8 m reflector, photoelectric photometer; MOST – the Canadian MOST satellite.

<table><caption>**Table 4.** Journal of the spectro-interferometric observations.</caption><thead><tr><th>Instr.</th><th>ΔT<br>(RJD)</th><th>ΔB<br>(m)</th><th>Δλ<br>(nm)</th><th>N<sub>V2</sub>/N<sub>T3φ</sub></th></tr></thead><tbody><tr><td>1</td><td>48 275–48 563</td><td>14–30</td><td>500–800</td><td>108/0</td></tr><tr><td>2</td><td>51 093–56 298</td><td>0–79</td><td>550–850</td><td>13 461/4 137</td></tr><tr><td>3</td><td>55 825–56 228</td><td>31–279</td><td>532–760</td><td>6 132/0</td></tr><tr><td>4</td><td>56 264–56 264</td><td>41–139</td><td>1 200–2 600</td><td>2 160/720</td></tr></tbody></table>

**Notes.** Δ*T* is the time span between the first and the last observation, Δ*B* the range of the projected baselines, Δ*λ* the wavelength range, *N*<sub>V2</sub> the total number of visibility observations, and *N*<sub>T3φ</sub> the total number of closure phase observations. In column “Instr.”: 1 – Mark III; 2 – NPOI; 3 – CHARA/VEGA; 4 – VLTI/AMBER.

The majority of the spectra at our disposal was acquired in three wavelength regions $Δλ \in \{4200–4500; 4750–5000; 6200–6700\}$ Å. Each region contains a Balmer line, which turned out to be the best for measuring the RVs of component B and several metallic lines, which gave accurate RVs of components Aa and Ab. These regions were also extracted from echelle spectra, and RVs were measured on each region independently. The last region ($Hα$) contains a number of telluric lines, including the $Hα$ line itself. Our model is unable to account for a telluric spectrum, and consequently it was not possible to measure accurate RVs of $Hα$ with this technique.

Initial RVs for the searching program were computed from the orbital solution presented in Nemravová et al. (2013), and we searched for the RV for each component in the interval $[-70; 70] \text{ km s}^{-1}$ that surrounds the initial estimate. The components of the eclipsing binary Aa and Ab are very similar, therefore we had to verify that the two components had not been interchanged by the program, especially near the conjunctions. If they were, the search was repeated using a narrower search interval.

The RVs and their uncertainty were estimated in the following way:

1. The parameters of synthetic spectra were chosen randomly from the Gaussian distributions centred at values listed in Table 7, and the standard deviations were set to their uncertainties.

2. The synthetic spectra were fitted to the observed ones. The procedure was repeated five hundred times for each spectrum, and the RV including its uncertainty was estimated from the resulting distribution.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Coverage of orbits 1 and 2 with the spectro-interferometric observations. The *outer plot*: the black line denotes the orbit of the centre of mass of the eclipsing binary relative to component B (which resides at the beginning of the coordinate system of the *outer plot*), and red dots denote the relative position of the centre of mass of the eclipsing binary relative to component B at the epochs of spectro-interferometric observations. The *inset plot*: the black line denotes the orbit of component Ab relative to component Aa (which resides at the beginning of the coordinate system of the *inset plot*), and red dots denote the relative position of component Ab relative to component Aa at epochs of spectro-interferometric observations. In both plots the orbital elements are invariable, i.e. they do not show the true orbits 1 and 2 as they would appear on the sky, but only demonstrate that the spectro-interferometric observations sample the orbits well enough to constrain elements of both orbits.</figcaption></figure>

This approach allowed us to estimate only the statistical part of the total uncertainty. The statistical uncertainty $ΔRV_{stat}$ was typically $\le 1 \text{ km s}^{-1}$ for components Aa and Ab and $\le 10 \text{ km s}^{-1}$ for component B. Measuring the RVs of component B was more difficult because the majority of metallic lines in its spectrum is very shallow and smeared out by the high rotational velocity of component B. The measurements are also very sensitive to the choice of the model and its discrepancies.

The telluric lines in the red and IR parts of the spectra were used to correct for the variations of the zero-point of the RV scale. These corrections were typically $\le 2 \text{ km s}^{-1}$ for the Ondřejov spectra, hence all measurements for which the RV
zero-point could not be checked in this way were assigned an uncertainty max(ΔRV<sub>stat</sub>, 2) km s<sup>-1</sup>, and the remaining ones were assigned an uncertainty max(ΔRV<sub>stat</sub>, 1) km s<sup>-1</sup>, where 1 km s<sup>-1</sup> is the upper bound of the precision of the zero-point correction for the Ondřejov spectra.

## 3.2. Direct analysis of RVs

Since we were not aware of any publicly available program for orbital solutions of hierarchical systems with apsidal advance of the outer orbit(s), JN developed such a program. The measured RVs were fitted with a model, which takes into account the two dynamical interactions between the three or four components. The effects considered are the apsidal motion of orbit 2 and the light-time (LITE) effect produced by orbits 2 (Δt<sub>LITE</sub> ≤ 0.006 d) and 3 (Δt<sub>LITE</sub> ≤ 0.013 d). The RVs of the jth component RV<sub>j</sub> were fitted with the standard Keplerian model:

$$RV_j(t) = \sum_i K_i [\cos(\omega_i(t) + v_i(t)) + e_i \cos \omega_i(t)], \quad (2)$$

where the index *i* goes over those orbits of ξ Tau that are relevant for the motion of the jth component of the ξ Tau system, $K_i$ is the semiamplitude of the RV curve, $\omega_i$ the argument of periastron, $v_i$ the true anomaly, $e_i$ the eccentricity, and *t* is time. The LITE $\Delta t_{LITE}$ was computed using the following formulae:

$$\Delta t_{\text{LITE},j}(t) = \sum_i \frac{P_i K_i (1 - e_i^2)^{\frac{3}{2}}}{2\pi c} \frac{\sin[\omega_i(t) + v_i(t)]}{1 + e_i \cos v_i(t)}, \quad (3)$$

where the index *i* goes over those orbits that are hierarchically above the orbit in which the jth component lies (i.e. over those that produce LITE), *P* is the orbital period, and *c* is the speed of light. Otherwise the notation is the same as for Eq. (2). The argument of periastron is a linear function of time $\omega_i(t) = \omega_i(t_0) + \dot{\omega}_i(t - t_{0,i})$, where $t_{0,i}$ is the reference epoch and $\dot{\omega}_i$ is the mean apsidal motion of the ith orbit.

The model elements were optimised by searching the minimum of the following $\chi^2$:

$$\chi^2 = \sum_{k=1}^{N_S} \sum_{j=1}^{N_C} \sum_{l=1}^{N_O} \frac{1}{\sigma_{j,l}^2} \left[ RV_j^{\text{OBS}}(\tilde{t}_{j,l}) - RV_j^{\text{SYN}}(\tilde{t}_{j,l}) - \gamma_k \right]^2, \quad (4)$$

where the index *k* goes over $N_S$ subsets of the measured RVs, which are defined in Table 2, the index *j* over $N_C$ components of the ξ Tau system for which RVs were measured, and the index *l* goes over $N_O$ individual measurements of the RV and $\tilde{t}$ is time corrected for the LITE. $\sigma$ denotes individual rms of the RVs estimated with the procedure described in Sect. 2, $RV^{\text{OBS}}$ the measured RV, $RV^{\text{SYN}}$ the model RV computed with Eq. (2), and corrected for the LITE via Eq. (3), and $\gamma$ denotes the systemic velocity. The minimum of the $\chi^2$ given by Eq. (4) was searched for with the sequential least-squares routine (Kraft 1988).

As discussed above, RVs of component B are less accurate than those of components Aa and Ab. Hence only RVs of the members of the eclipsing binary were fitted to obtain the majority of orbital elements. The individual subsets for individual types of the spectra gave very similar values of the systemic velocity (within 3$\sigma$), hence all available measurements were grouped together and a joint systemic velocity was derived for them. When a final solution was obtained, the measurements were complemented with RV measurements of component B and the mass ratio $q_2$ was optimised (keeping the remaining parameters fixed). The parameters corresponding to the best-fit solution

<table><caption>**Table 5.** Parameters of the two-orbit (1 and 2) fit given by Eqs. (2) and (3) to measured RVs.</caption><thead><tr><th rowspan="2">El. Orb.</th><th rowspan="2">Units</th><th colspan="2">Values</th></tr><tr><th>1</th><th>2</th></tr></thead><tbody><tr><td>P<sub>AN</sub></td><td>(d)</td><td>—</td><td>145.579 ± 0.048</td></tr><tr><td>P<sub>S</sub></td><td>(d)</td><td>7.14664 ± 0.00002</td><td>145.113 ± 0.071</td></tr><tr><td>T<sub>min</sub></td><td>(RJD-56 220)</td><td>4.7067 ± 0.0025</td><td>—</td></tr><tr><td>T<sub>p</sub></td><td>(RJD-55 600)</td><td>—</td><td>9.46 ± 0.52</td></tr><tr><td>K</td><td>(km s<sup>-1</sup>)</td><td>87.79 ± 0.25</td><td>38.37 ± 0.19</td></tr><tr><td>e</td><td></td><td>0.0<sup>1</sup></td><td>0.2101 ± 0.0053</td></tr><tr><td>q</td><td></td><td>0.9438 ± 0.0036</td><td>0.889 ± 0.056</td></tr><tr><td>ω</td><td>(deg)</td><td>90<sup>1</sup></td><td>9.25 ± 1.42</td></tr><tr><td>ω̇</td><td>(deg yr<sup>-1</sup>)</td><td>0.0<sup>1</sup></td><td>2.90 ± 0.33</td></tr><tr><td>N</td><td></td><td>748</td><td></td></tr><tr><td>χ<sup>2</sup><sub>R</sub></td><td></td><td>2.128</td><td></td></tr><tr><td colspan="4" style="text-align: center;">Systemic velocity</td></tr><tr><td>γ</td><td>(km s<sup>-1</sup>)</td><td colspan="2">8.05 ± 0.18</td></tr></tbody></table>

**Notes.** $P_{AN}$ denotes the anomalous period, $P_S$ the sidereal period, $T_{min}$ the epoch of the primary minimum of the light curve, and $T_p$ the epoch of the periastron passage. (1) The parameter was fixed. $K_1$ refers to the primary of the eclipsing binary $K^{Aa}$, and $K_2$ to the centre of gravity of the eclipsing binary $K^{Aa+Ab}$.

are listed in Table 5. RVs and the best-fitting model are plotted against time (to show the secular evolution of the periastron argument) for orbit 2 in Fig. 2, and against phase for orbit 1 in Fig. 3.

The uncertainties and correlations of individual parameters were estimated with the bootstrap method. One thousand samples were randomly chosen from all available RVs. Each sample consisted of the same number (748) of measurements as the original (meaning that some measurements repeat within a sample). Each sample was fitted with an orbital model and the uncertainties were estimated from the distribution of the results.

The reduced $\chi^2$ (denoted $\chi_R^2$ throughout the article) $\chi_R^2 \approx 2$, which is greater than ideal case of 1, is probably caused by variations of the RV zero-point larger than we accounted for (we note that the estimate is based on the variations of the zero-point measured on the Ondřejov red spectra), and by the fact that the synthetic spectra need not correspond to the observed ones in all details, for which we cannot account properly. Moreover, the model does not account properly for the dynamical interaction (see Sects. 8 and 9) between all orbits.

We also fitted a model including orbit C fixed at the orbital elements given in Table 10. The reduced $\chi^2$ was only marginally (≤1%) lower than that in Table 5. This is expected because the semi-amplitude of the RV caused by the revolution of the triple subsystem around the common centre of gravity with component C is ≈1 km s<sup>-1</sup> and the LITE produced by that motion is ≈0.013 d, which means that both are beyond the detection limit of our measurements.

## 3.3. Spectral disentangling

We were only able to separate the spectra in the vicinity of five major spectral lines Hα, Hβ, He I 4471 Å, Mg II 4481 Å, and Hγ because only these regions were available for both the slit and echelle spectra. An attempt was made to separate the spectra of individual components using only the spectra from the three available echelle spectrographs. However, these separated spectra had strongly warped continua and were unsuitable for further investigation.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** RVs of the centre of gravity of the eclipsing binary (red triangles) and component B (blue triangles) against the best-fitting model (black) corresponding to parameters listed in Table 5. $\Delta_{Aa,Ab}$ (in km s$^{-1}$) denote residuals of the fit for RVs of the centre of gravity of the eclipsing binary, and $\Delta_B$ (in km s$^{-1}$) residuals of the fit for RVs of component B.</figcaption></figure>

We used the program KOREL (Hadrava 1995, 1997, 2009) (release 04-2004), which not only separates the spectra, but also fits the spectroscopic orbital elements. This gave us the opportunity to compare the orbital solution obtained directly from the measured RVs with the result of KOREL. Only components B, Aa, and Ab were fitted because component C is not detectable. Relative luminosities of all three components were kept constant during the orbital motion. This assumption, although not exactly satisfied because of the shallow eclipses of components Aa and Ab, was necessary for the stability of the disentangling.

The orbital elements presented in Table 5 served as the starting estimates for the minimisation. The spectroscopic orbital elements obtained with KOREL are listed in Table 6. The separated profiles from the considered spectral regions are shown in Fig. A.1. KOREL does not provide the uncertainties of the fitted elements. Therefore a map of the $\chi^2$ around the minimum found with the minimisation engine was drawn for every combination of two fitted parameters. The uncertainties, which are listed in Table 6, correspond to 68% confidence intervals (roughly one $\sigma$) estimated from these maps.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** RVs of components Aa (red) and Ab (blue) relative to the centre of gravity of the eclipsing binary against the best-fitting model (black) listed in Table 5. $\Delta_{Aa,Ab}$ are residuals of the fit for components Aa and Ab.</figcaption></figure>

<table><caption>**Table 6.** Orbital elements obtained by KOREL (spectral disentangling) for all available spectra containing at least one of the studied regions.</caption><thead><tr><th>Elem. orbit</th><th>Unit</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>P<sub>AN</sub></td><td>(d)</td><td>—</td><td>145.612 ± 0.056</td></tr><tr><td>P<sub>S</sub></td><td>(d)</td><td>7.14664 ± 0.00002</td><td>145.123 ± 0.072</td></tr><tr><td>T<sub>min</sub></td><td>(RJD-56 220)</td><td>4.6963 ± 0.0040</td><td>—</td></tr><tr><td>T<sub>p</sub></td><td>(RJD-56 000)</td><td>—</td><td>9.29 ± 1.44</td></tr><tr><td>K</td><td>(km s<sup>-1</sup>)</td><td>87.52 ± 0.59</td><td>37.55 ± 0.57</td></tr><tr><td>e</td><td></td><td>0<sup>1</sup></td><td>0.180 ± 0.024</td></tr><tr><td>q</td><td></td><td>0.943 ± 0.008</td><td>1.02 ± 0.27</td></tr><tr><td>ω</td><td>(deg)</td><td>90<sup>1</sup></td><td>8.52 ± 4.1</td></tr><tr><td>ω̇</td><td>(deg yr<sup>-1</sup>)</td><td>0<sup>1</sup></td><td>3.032 ± 0.38</td></tr><tr><td>χ<sub>R</sub><sup>2</sup></td><td></td><td></td><td>1.19</td></tr></tbody></table>

**Notes.** The orbital model consists of orbits 1 and 2. (1) The parameter was fixed. $K_1$ refers to the primary of the eclipsing binary $K^{Aa}$ and $K_2$ to the centre of gravity of the eclipsing binary $K^{Aa+Ab}$.

An attempt was carried out to separate the lines of component C within two spectral bands in the near infrared, $\Delta\lambda_{IR} = \{7750-7800, 8570-8800\} \AA$. The spectrum of component C was not detected in either of these bands. It was probably caused by the relatively low signal-to-noise ratio (S/N) of the echelle spectra in the infrared region and their limited number.

We also note that we tried to use the separated profiles instead of synthetic ones to measure RVs with the PYTERPOL program written by JN. This worked well for components Aa and Ab, but failed for component B. The reason is that the shape of the separated spectral lines depends on the orbital elements, for which the spectra were separated, and vice versa. Hence the separated spectra partially “remember” the orbital elements for which they were obtained, and if they are used for the RV measurements, they would give a fine RV curve described by a solution close to these elements. This becomes a problem when one or more orbital elements suffer from a large uncertainty, which was the case for $\xi$ Tau in the mass ratio of orbit 2.

### 3.4. Comparison of observed and synthetic spectra

JN has developed a Python program PYTERPOL<fnref n="2" />, which interpolates in a pre-calculated grid of synthetic spectra to obtain estimates of the radiative properties of the components of multiple systems. For $\xi$ Tau these parameters were the effective temperature $T_{\text{eff}}$, gravitational acceleration $\log g$, the projected rotational velocity $v \sin i$, RV, and the relative luminosity $L_R$. The parameters of components Aa, Ab, and B were covered by the POLLUX grid (Palacios et al. 2010), and component C was searched for using the AMBRE grid (de Laverny et al. 2012). Solar metallicity was assumed.

The fit was carried out in four spectral regions, but only three relative luminosities were derived, since two of the regions are very close to each other and the luminosities $L_R$ are most likely almost the same.

The spectral regions were $\Delta\lambda_1 = [4280, 4495] \AA$, $\Delta\lambda_2 = [4815, 4940] \AA$, and $\Delta\lambda_3 = \{[6330, 6390]; [6660, 6695]\} \AA$.

The relative luminosities were assumed to be constant over each spectral region $\Delta\lambda_i$.

Two of the regions contain a Balmer line, which constrains the gravitational acceleration of all three components, and also a large number of metallic lines, which constrain the temperature, RVs, and the projected rotational velocities. We fitted 137 spectra from the Ondřejov Observatory together because their normalisation is straightforward (a first-order polynomial often suffices to fit the continuum), so that the Balmer lines are not affected by systematics often introduced by the normalisation. The uncertainty of the relative flux was estimated from the continuum for each spectrum and set constant for each spectrum.

The bootstrap method was used to obtain a best-fit set of parameters. We randomly drew 137 spectra from the pool of 137 Ondřejov spectra (meaning that one or more spectra can be present multiple times within the random sample) and fitted them. The initial set of parameters was randomly chosen from intervals<fnref n="3" /> which were established from the first trial fits. The initial RVs were estimated from the orbital solution presented in Nemravová et al. (2013) and randomly put slightly off (within 30 km s⁻¹ vicinity of the estimate) to secure robustness of the final solution. The procedure was repeated five hundred times and the final set of parameters was estimated from the distribution of the results. The shape of the distribution was Gaussian-like, that is, describable with a mean value and its standard deviation. The results are presented in Table 7.

A comparison of four spectral regions with the model is shown in Fig. 4. The reduced $\chi_R^2$ is lower than one, indicating that we have slightly overestimated the uncertainty of the relative flux of the observed spectra.

### 3.5. Comparison of synthetic and separated spectra

We fitted the separated spectra corresponding to the solution of Table 6 with the interpolated synthetic spectra to check the results of Sect. 3.4. The program PYTERPOL was used again.

<footnote n="2">² A detailed description with a simple tutorial how to use it is provided at https://github.com/chrysante87/pyterpol/wiki</footnote>

<footnote n="3">³ The intervals are the following: $T_{\text{eff}}^{\text{B}} \in [13\,000, 14\,500]$ K, $T_{\text{eff}}^{\text{Aa}} \in [9000, 11\,500]$ K, $T_{\text{eff}}^{\text{Ab}} \in [9000, 11\,500]$ K, $\log g^{\text{B}} \in [4.0, 5.0]$, $\log g^{\text{Aa}} \in [3.5, 4.5]$, $\log g^{\text{Ab}} \in [3.5, 4.5]$, $v \sin i^{\text{B}} \in [200, 250]$ km s⁻¹, $v \sin i^{\text{Aa}} \in [0, 40]$ km s⁻¹, $v \sin i^{\text{Ab}} \in [0, 40]$ km s⁻¹, $L_R^{\text{B}} \in [0.55, 0.8]$, $L_R^{\text{Aa}} \in [0.10, 0.25]$, $L_R^{\text{Ab}} \in [0.10, 0.25]$.</footnote>
<table><caption>**Table 7.** Parameters of the fit of the synthetic spectra to 137 observed Ondřejov spectra.</caption><thead><tr><th>Parameter<br>Component</th><th>Unit</th><th>B</th><th>Aa</th><th>Ab</th></tr></thead><tbody><tr><td>T<sub>eff</sub></td><td>(K)</td><td>14 190 ± 150</td><td>10 700 ± 160</td><td>10 480 ± 130</td></tr><tr><td>log g</td><td>(cgs)</td><td>4.527 ± 0.041</td><td>4.08 ± 0.12</td><td>4.01 ± 0.10</td></tr><tr><td>v sin i</td><td>(km s<sup>-1</sup>)</td><td>229.2 ± 1.7</td><td>12.6 ± 2.6</td><td>14.3 ± 3.1</td></tr><tr><td>L<sub>R</sub><sup>Δλ<sub>1</sub></sup></td><td></td><td>0.660 ± 0.024</td><td>0.179 ± 0.018</td><td>0.165 ± 0.022</td></tr><tr><td>L<sub>R</sub><sup>Δλ<sub>2</sub></sup></td><td></td><td>0.688 ± 0.026</td><td>0.162 ± 0.024</td><td>0.155 ± 0.027</td></tr><tr><td>L<sub>R</sub><sup>Δλ<sub>3</sub></sup></td><td></td><td>0.665 ± 0.036</td><td>0.173 ± 0.028</td><td>0.161 ± 0.031</td></tr><tr><td>χ<sub>R</sub><sup>2</sup></td><td></td><td></td><td>0.87</td><td></td></tr></tbody></table>

**Notes.** The modelled spectral intervals are $\Delta\lambda_1 = [4280, 4495]$ Å, $\Delta\lambda_2 = [4815, 4940]$ Å and $\Delta\lambda_3 = {[6330, 6390]; [6660, 6695]}$ Å.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Example of the fit of the synthetic spectra (red) to three observed spectra (black) in spectral regions: 1) $\Delta\lambda_1 = [4280, 4495]$ Å (top), 2) $\Delta\lambda_2 = [4815, 4940]$ Å (middle), 3) $\Delta\lambda_3 = [6330, 6390]$ Å (bottom, left), 4) $\Delta\lambda_3 = [6660, 6695]$ Å (bottom, right). The synthetic spectra are given by parameters listed in Table 7.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Fit of the light curve from the satellite MOST. Only the light curve minima and their surroundings are shown. The primary (secondary) minimum is on the left (right) on each panel. The left panel corresponds to the global circular solution $e_1 = 0.0$ and to orbital period $P_1 = 7.14664$ d. The right panel corresponds to a local solution, where small adjustment of the eccentricity and the orbital period was allowed. MO denotes the satellite broad-band filter.</figcaption></figure>

<table><caption>**Table 8.** Parameters of synthetic spectra best-fitting the separated spectra.</caption><thead><tr><th rowspan="2">Element</th><th rowspan="2">Unit</th><th colspan="3">Value component</th></tr><tr><th>B</th><th>Aa</th><th>Ab</th></tr></thead><tbody><tr><td>$T_{\text{eff}}$</td><td>(kK)</td><td>14.07(14)</td><td>10.26(14)</td><td>10.050(80)</td></tr><tr><td>$\log g_{[\text{cgs}]}$</td><td></td><td>3.99(4)</td><td>4.06(9)</td><td>4.02(4)</td></tr><tr><td>$v \sin i$</td><td>(km s<sup>-1</sup>)</td><td>253.6(16)</td><td>18.6(12)</td><td>10.2(10)</td></tr><tr><td>$L_{\text{R}}^{\Delta\lambda_1}$</td><td></td><td>0.758(8)</td><td>0.168(3)</td><td>0.150(7)</td></tr><tr><td>$L_{\text{R}}^{\Delta\lambda_2}$</td><td></td><td>0.711(6)</td><td>0.191(5)</td><td>0.149(3)</td></tr><tr><td>$L_{\text{R}}^{\Delta\lambda_3}$</td><td></td><td>0.686(7)</td><td>0.188(4)</td><td>0.161(7)</td></tr><tr><td>$\gamma$</td><td>(km s<sup>-1</sup>)</td><td></td><td>8.1(27)</td><td></td></tr><tr><td>$\chi_{\text{R}}^2$</td><td></td><td></td><td>31.58</td><td></td></tr></tbody></table>

**Notes.** $\gamma$ denotes the systemic velocity of $\xi$ Tau. The fit is plotted in Fig. A.1.

The following spectral regions were fitted:

$$ \Delta\lambda_1 = \{[4280, 4400]; [4455, 4495]\} \AA, $$

$$ \Delta\lambda_2 = [4765, 4970] \AA, \text{ and} $$

$$ \Delta\lambda_3 = \{[6325, 6395]; [6510, 6620]; [6655, 6695]\} \AA. $$

The parameters corresponding to the best-fitting synthetic spectra are listed in Table 8. The best-fit parameters were estimated with a MCMC simulation and the uncertainties reflect only the statistical part of the uncertainty. The systematic uncertainty – the warp in the continua and the need for its normalisation – cannot be easily quantified and is responsible for the extremely high reduced $\chi^2$ along with the very high S/N ratio of the separated spectra. Therefore the uncertainties of the parameters listed in Table 8 are very likely underestimated.

This systematic effect corrupts the estimate of $\log g$ of all components, especially component B, where the warping was the most pronounced; therefore it also applies to the rotational velocity of component B. The rotational velocity of

components Aa and Ab is strongly affected by the choice of the instrumental broadening, which is very difficult to estimate for separated spectra and was set to 0.2 Å. The total light is also very likely affected by the re-normalisation, which (necessarily) changes the depths of spectral lines ($L = \sum_{i=1}^{3} L_{\text{R}}^{\Delta\lambda_{i}} \neq 1$ for all studied bands).

Bearing all this in mind, we state that this result does not contradict, but rather supports that obtained by fitting of synthetic to observed spectra. A comparison of the synthetic spectra corresponding to the parameters listed in Table 8, of separated spectra, and of re-normalised separated spectra is in Fig. A.1.

## 4. Photometry

The preliminary analysis published in Nemravová et al. (2013) has shown that the light variations can be attributed to the eclipses of components Aa and Ab of orbit 1. They partially eclipse each other and produce two very narrow and nearly identical minima, which are only ≈0.1 mag deep in the Johnson V passband.

In addition to the binary eclipses, our new very precise MOST satellite observations unveiled persistent low-amplitude rapid cyclic light changes that are probably associated with component B, since they remain during both binary eclipses. The MOST light curve also allows determining very accurate radii of components Aa and Ab as well as detecting variations of the mean motion of the eclipsing pair. The zoomed parts of both minima of the MOST light curve are shown in Fig. 5.

### 4.1. Period analysis of the light curve

Our first goal in the analysis of the MOST light curve was to unveil the nature of the rapid cyclic low-amplitude changes. Two different methods were used to construct and investigate the periodogram of the light curve. The first is based on the Fourier transform (FT hereafter) and is implemented in the program PERIOD04 (Lenz & Breger 2004). The second uses the
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Fourier spectrum of the MOST light curve from Fig. 7 (i.e. out-side eclipses). Prominent frequency peaks are marked (see their descrip-tion in Sect. 4.1).</figcaption></figure>

phase dispersion minimisation technique (PDM) (Stellingwerf 1978) and is implemented in the program HEC27⁴. The peri-odogram of the whole light curve is dominated by the orbital period of the eclipsing binary $P_1 \approx 7.147$ d. To study the rapid low-amplitude oscillation, we removed the eclipses (see Fig. 7, top).

The periodogram of the rapid oscillations (see Fig. 6) shows a basic frequency of $f_0 = 2.38$ d⁻¹, most likely due to rotation of component B, the first harmonics of the eclipsing binary or-bital frequency $f_1 = 2/P_1 = 0.279$ d⁻¹; the frequencies of $f_d = 1.002738$ d⁻¹ and $f_{MOSTorbit} = 14.2$ d⁻¹ are instrumental (i.e. the orbital frequency of the satellite). The remaining prominent fre-quencies $f_{alias} = \{15.1734, 17.5385, 28.3896, 42.5825, 56.7745, 70.9720\}$ d⁻¹ seem to be either integer multiples of $f_{orb}$ or its splittings with $f_0$ or $f_d$. Remaining peaks (e.g. $f = 87.1609$ d⁻¹) have relatively low S/N ratios. We are not aware of any instru-mental effect that would induce oscillations at $f_0 = 2.38$ d⁻¹, hence the low-amplitude variations arise from a physical process in $\xi$ Tau.

A closer look at Fig. 7 shows that the amplitude of the curve varies. To quantify these changes, a harmonic function $f(t) = 1 + C_0 + A_0 \sin[2\pi(t - T_0)f_0 + \phi_0]$ was sequentially fitted to segments of the light curve $\Delta t_1 = P_1/2$ d wide, and shifted with a step $\Delta t_2 = P_1/20$, where $P_1$ is the period of the eclipsing binary. The scan revealed that both the basic frequency $f_0$ and its amplitude $A_0$ vary on the time span of two orbital periods of the eclipsing binary (see Fig. 7, middle and bottom panels).

## 4.2. Nature of quasiperiodic oscillations

The quasiperiodic oscillations clearly visible in the MOST light curve with an approximate period $P_0 \approx (0.42 \pm 0.01)$ d and an amplitude $A_0 = (0.00060 \pm 0.00015)$ mag exhibit both a fre-quency (FM) and an amplitude modulation (AM) on the time span of about the two shortest orbital periods $P_1$ (see Fig. 7). We can think of several possibilities regarding their origin: an instru-mental effect, a fifth component and ellipsoidal variations, rota-tion with spots, or rotation and pulsations.

The first option does not seem very likely, however, because we do not know about any instrumental period of 0.42 d (like one day, or a satellite orbital period 0.07042 d in this case).

A hypothetical fifth component (second option) orbiting ei-ther component B, Aa, or Ab with a period $2P_0$ can induce ellipsoidal variations of the order of $A_0$, but they would be expected to be very regular (without large AM, FM) and to manifest themselves in one of the RV curves as well, which is not the case. We do not see any peak in the Fourier spectrum at $f_0 = 1/P_0 = 2.38$ d⁻¹, even though the Nyquist frequency for our spectroscopic dataset is $f_{Ny} = 7.1$ d⁻¹. Nevertheless, the coverage and cadence are not uniform at all and the expected am-plitude is small (5 km s⁻¹), which makes this particular argument weak. We would also expect to see some frequency modulation due to the (classical) Doppler effect, $P'_{obs} = (1 - \frac{v}{c})\gamma P_{src}$, with $v \approx 2v_{kepl}$. However, for 0.423 d we would only obtain a change by 0.001 d, which is one order of magnitude smaller than the observed total variation.

The lower limit for the rotation period is the critical rota-tion, $P_{min} = 2\pi(GM/R^3)^{-1/2}$, and the upper limit is determined by rotational broadening, $P_{max} = 2\pi R/(v \sin i)$ (cf. Table 8). For component Aa or Ab, the admissible range is from about $P_{rot} = 0.180$ d to 3.85 d, for component B it is 0.325 to 0.634 d. The observed oscillations are within both ranges, so that we can-not distinguish the source component at this point. One can argue that small axial inclination for components Aa, Ab is unlikely when their orbital inclination is large, so that their true $P_{rot} > P_0$. We thus prefer to attribute these oscillations to component B. Additionally, this star is relatively brighter so that it is easier to induce the oscillations of given amplitude $A_0$.

It seems difficult to distinguish between spots and pulsations (options three and four above; as in Degroote et al. 2011). Especially for early-type stars, spots are infrequent, unless a star is chemically peculiar or magnetically active (Bp), but we have no observations and analyses at our disposal that could prove or disprove this for $\xi$ Tau.

Pulsating B stars (like $\beta$ Cep, SPB) always exhibit a low-frequency signal corresponding to the rotation and then a series of pulsation modes, either pressure (high-frequency) or grav-ity (low-frequency). The cadence of MOST photometric ob-servations allows us to compute the Fourier spectrum up to $f_{Ny} = 719$ d⁻¹, corresponding to 0.00139 d = 2 min (Fig. 6). Except for the basic rotational period, its aliases with the orbital period $P_1$ of the eclipsing binary, one-day and $P_{orb}$ instrumen-tal periods, we can unfortunately not unambiguously detect any pulsation modes with $S/N \ge 5$, to say nothing about rotational splittings, which would be conclusive.

## 4.3. Eclipse timing variations

The orbital period of the eclipsing binary $P_1 = 7.14664$ d in-troduces a small but clearly detectable shift $\Delta_{PHASE} \approx 0.0003$ between the two minima recorded with the MOST satellite. The shift disappears if the orbital period and the eccentricity are opti-mised. The local period and eccentricity, which do not cause the phase shift, are $P_1 = 7.14466$ d and $e_1 \approx 0.002$. The problem is illustrated in Fig. 5, where the comparison of an eccentric model with the local value of the orbital period and a global circular model is shown. An even larger phase shift $\Delta p \sim 0.004$ was de-tected when a similar analysis was carried out for all photometric observations.

This led us to investigate the eclipse timing variations (ETVs) in all available photometry, divided into subsets covering time intervals shorter than $P_2/4$ (individual minima are shown in Figs. B.1 and B.2). The ETVs are very noisy, and the delays themselves have an amplitude $\Delta t_{OBS} \approx 0.025 \pm 0.01$ d that cannot be explained by LITE ($\Delta t_{LITE} \le 0.006$ d). Moreover, they seem

<footnote n="4">⁴ The program and a short user guide are available at http://astro.troja.mff.cuni.cz/ftp/hec/HEC27</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Normalised light curve as reduced from MOST photometry, but without intervals of primary and secondary eclipses (top panel), together with the corresponding period $P_0$ (middle) and amplitude $A_0$ (bottom) of the harmonic function $f(t) = 1 + C_0 + A_0 \sin[2\pi(t - T_0)/P_0 + \phi_0]$, which was sequentially fitted to the light curve, always in limited intervals $\Delta E_1 = 0.5$ of the epoch (indicated by the black double arrow), shifted with a step $\Delta E_2 = 0.05$. The oscillations exhibit both frequency and amplitude modulations, with periods spanning $P_0 = (0.42 \pm 0.01)$ d and amplitudes $A_0 = (0.00060 \pm 0.00015)$ mag. It seems that the longest $P_0$ and the largest $A_0$ are observed at around primary eclipses and vice versa.</figcaption></figure>

to vary on a timescale comparable to the orbital period $P_2$. Hence we assume that the dynamical interaction between orbits 1 and 2 is the reason for these delays. The first-order model of the physical delay (Eq. (8) from Rappaport et al. 2013), which is only a part of the total ETV, arising from dynamical interaction of two orbits in hierarchical triple systems, gives an estimate of the amplitude of the effect $\Delta t_{\text{MODEL}} \approx 0.02$ d, (i.e. in rough agreement with the detected value). This is another proof of the dynamical interaction in ξ Tau (the first is the apsidal motion reported by Nemravová et al. 2013) and led us to develop an N-body model (see Sect. 8) and a perturbation theory (see Sect. 9).

## 4.4. Global orbital model for all light curves

The program PHOEBE 1.0 (Prša & Zwitter 2005, 2006) was used to derive the light-curve solution. The mass ratio $q_1$ was taken from the analysis of the RVs (see Table 5) because only light curves were modelled and they do not constrain the mass ratio for a detached system. The eccentricity was assumed to be $e_1 = 0.0$ (although Sect. 8 shows that orbit 1 is slightly eccentric). The value of the semi-major axis $a$ was adjusted after each iteration based on $a_1 \sin i$ given by the fit of the directly measured RVs (see Table 5). The linear limb-darkening law was adopted and the coefficients were interpolated in a pre-calculated grid distributed along with PHOEBE. The bolometric albedos were taken from Claret (2001) and the gravity brightening coefficients from Claret (1998) for the corresponding temperatures of components of the eclipsing binary. The spin-orbit synchronisation, that is, the synchronicity ratios $F^{\text{Aa}} = F^{\text{Ab}} = 1$, was assumed, because radii $R^{\text{Aa}}$ and $R^{\text{Ab}}$ from Nemravová et al. (2013) and rotational velocities from Table 7 give synchronicity ratios $F^{\text{Aa}} = 1.12 \pm 0.26$, and $F^{\text{Ab}} = 0.74 \pm 0.20$; the deviations from the corotation are small and probably arise from an incorrect determination of the radii. The primary effective

temperature $T_{\text{eff}}^{\text{Aa}}$ was set to the value found through a comparison of synthetic and observed spectra.

The orbital inclination $i_1$, Kopal surface potentials $\Omega_K^{\text{Aa}}$, $\Omega_K^{\text{Ab}}$ of both components and the epoch of the primary minimum $T_{\text{min},1}$, the secondary temperature $T_{\text{eff}}^{\text{Aa}}$, and the relative luminosity of component B $L^B$ in each spectral band were optimised. Initial estimates of these parameters were taken from Nemravová et al. (2013), initial relative luminosities $L^B$ of component B were estimated from the comparison of synthetic and observed profiles (Table 7). The primary luminosities $L^{\text{Ab}}$ were adjusted after each iteration.

The fitting was carried out in the Python environment of PHOEBE, and the minimum was determined with the differential evolution algorithm (Storn & Price 1997). The following parametric space was searched: $T_{\text{min},1} \in [56\,224.68, 56\,224.78]$ RJD, $i_1 \in [84, 90]$ deg, $\Omega_K^{\text{Aa}} \in [11, 20]$, $\Omega_K^{\text{Ab}} \in [11, 20]$, $T_{\text{eff}}^{\text{Ab}} \in [10\,000, 10\,700]$ K, $L^B \in [0.55, 0.78]$. The last interval applies to each studied spectral filter (U, B, V, MOST). The parametric space was densely sampled with models during the fitting ($\approx 300\,000$ light curve models were computed). This showed that the relative luminosity of component B $L^B$ is poorly constrained.

After a global minimum was found, we split our data and optimised the ephemeris, relative luminosity of component B, and surface potentials using only observations from the MOST satellite, after which we optimised the effective temperature of component Ab and the relative luminosity of component B using the Johnson UBV photometry. The epoch of the primary minimum was also fitted for the UBV dataset to slightly adapt it for the ETVs discussed in Sect. 4.3.

The parameters corresponding to the best-fitting model are listed in Table 9. Our model is unable to account for either the rapid light oscillations or the ETVs; therefore we raised the uncertainty of observations from the MOST satellite to deal with
<table><caption>**Table 9.** Parameters of the best-fitting circular orbital model obtained with the program PHOEBE 1.0.</caption><thead><tr><th>Element</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Orbital properties</td></tr><tr><td>P</td><td>(d)</td><td>7.14664 ± 0.00010</td></tr><tr><td>T<sub>min</sub></td><td>(RJD)</td><td>56 224.72482 ± 0.00022</td></tr><tr><td>a</td><td>(R<sub>☉</sub>)</td><td>25.552 ± 0.097</td></tr><tr><td>q</td><td></td><td>0.9439<sup>1</sup></td></tr><tr><td>e</td><td></td><td>0.0<sup>1</sup></td></tr><tr><td>i</td><td>(deg)</td><td>86.85 ± 0.22</td></tr><tr><td>ω</td><td>(deg)</td><td>90</td></tr><tr><td colspan="3">Component properties</td></tr><tr><td>Comp.</td><td></td><td>Aa</td><td>Ab</td></tr><tr><td>T<sub>eff</sub></td><td>(K)</td><td>10 700<sup>1</sup></td><td>10 450 ± 150</td></tr><tr><td>Ω<sub>K</sub></td><td></td><td>15.97 ± 0.25</td><td>15.93 ± 0.32</td></tr><tr><td>L<sub>V</sub></td><td></td><td>0.204 ± 0.020</td><td>0.174 ± 0.017</td></tr><tr><td>L<sub>B</sub></td><td></td><td>0.196 ± 0.021</td><td>0.165 ± 0.018</td></tr><tr><td>L<sub>U</sub></td><td></td><td>0.157 ± 0.012</td><td>0.130 ± 0.010</td></tr><tr><td>L<sub>MO</sub></td><td></td><td>0.203 ± 0.008</td><td>0.162 ± 0.006</td></tr><tr><td>L<sub>Hp</sub></td><td></td><td>0.212 ± 0.023</td><td>0.180 ± 0.020</td></tr><tr><td colspan="4">Passband luminosity of component B</td></tr><tr><td>L<sub>V</sub><sup>B</sup></td><td></td><td>0.622 ± 0.060</td><td></td></tr><tr><td>L<sub>B</sub><sup>B</sup></td><td></td><td>0.639 ± 0.069</td><td></td></tr><tr><td>L<sub>U</sub><sup>B</sup></td><td></td><td>0.713 ± 0.071</td><td></td></tr><tr><td>L<sub>MO</sub><sup>B</sup></td><td></td><td>0.634 ± 0.024</td><td></td></tr><tr><td>L<sub>Hp</sub><sup>B</sup></td><td></td><td>0.608 ± 0.067</td><td></td></tr><tr><td>χ<sub>R</sub><sup>2</sup></td><td></td><td>1.134</td><td></td></tr></tbody></table>

**Notes.** All available photometric observations were fitted. $\Omega_K$ denotes the Kopal surface potential here, and $L$ the relative luminosity in the filter given by the subscript. U, B, and V denote Johnson filters, MO denotes the broad-band filter in the MOST satellite, and Hp denotes the broad-band filter of the satellite HIPPARCOS. (1) The parameter was kept fixed.

the former ($\Delta m_{MOST} = 0.006$ given by the sinusoidal fit). The uncertainties of parameters are estimated as 68% confidence intervals computed from a scaled $\chi^2$ (scaled to an ideal situation, where the $\chi^2_R = 1$), although in this case the scaling was almost unnecessary, since the best solution has $\chi^2_R = 1.134$.

## 5. Astrometry of orbit 3

We used the existing astrometric positions listed in the WDS catalogue (see Mason et al. 1999, and references therein) to improve the orbital elements of orbit 3 published by Rica Romero (2010). The solution was carried out with the help of the program written by PZ (see Zasche & Wolf 2007, and references therein). The solution is listed in Table 10 and the orbit is shown in Fig. 8.

## 6. Spectro-interferometry

In this section we present an orbital analytic model of the $\xi$ Tau system, which we fit to spectro-interferometric observations to estimate orbital elements, radii, and fractional luminosities of $\xi$ Tau.

<table><caption>**Table 10.** Orbital elements of orbit 3 based on a fit to astrometric measurements published in WDS.</caption><thead><tr><th>Element</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>P</td><td>(yr)</td><td>51.01 ± 0.78</td></tr><tr><td>T<sub>p</sub></td><td>(RJD)</td><td>54 615 ± 251</td></tr><tr><td>e</td><td></td><td>0.5728 ± 0.0028</td></tr><tr><td>a</td><td>(mas)</td><td>441.5 ± 2.4</td></tr><tr><td>i</td><td>(deg)</td><td>25.4<sup>1</sup> ± 7.7</td></tr><tr><td>ω</td><td>(deg)</td><td>10.6 ± 8.9</td></tr><tr><td>Ω</td><td>(deg)</td><td>106.4 ± 2.2</td></tr></tbody></table>

**Notes.** The listed parameters are the orbital period $P$, the periastron epoch $T_p$, the eccentricity $e$, the semi-major axis $a$, the inclination $i$, the argument of the periastron $\omega$, and the position angle of the nodal line $\Omega$. (1) The inclination is determined ambiguously. A solution with $i = -25.4$ deg has exactly the same $\chi^2$.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Speckle-interferometric outer orbit 3 corresponding to the solution of Table 10. The dotted line stands for the line of apsides, the dashed line for the line of nodes.</figcaption></figure>

### 6.1. Global model for all available spectro-interferometric observations

The calibrated visibilities from VEGA/CHARA were fitted night by night with a model consisting of three uniform disks using the tool LitPro<fnref n="5" /> (Tallon-Bosc et al. 2008). The observations obtained during each single night were not numerous enough to safely estimate the positions and radii of components Aa, Ab, and B on the celestial sphere. In contrast to this, the NPOI observations are numerous enough to provide good estimates of the relative position of component B and the photocentre of the eclipsing binary for each night. They are presented in Table C.1 along with details on their acquisition (see Appendix C).

To circumvent the problem, we created a global orbital model that computes instantaneous positions of components B,

<footnote n="5">⁵ LITpro software available at http://www.jmmc.fr/litpro</footnote>
Aa, and Ab with the following formulae:

$$ \alpha_i(t) = \arctan(\tan(v_i(t) + \omega_i(t)) \cos i_i) + \Omega_i, \quad (5) $$

$$ \rho_i(t) = a_i \frac{1 - e_i^2}{1 + e_i \cos v_i(t)} \frac{\cos(\omega_i(t) + v_i(t))}{\cos(\alpha_i(t) - \Omega_i)}, \quad (6) $$

$$ x_i = \rho_i \sin \alpha_i, \quad (7) $$

$$ y_i = \rho_i \cos \alpha_i, \quad (8) $$

where index *i* denotes the component of a binary, *v* is the true anomaly, *ω* the argument of periastron, *i* the orbital inclination with respect to the celestial sphere, *Ω* is the position angle of the nodal line, *a* the angular semimajor axis, and *e* the eccentricity. The position angle *α*<sub>*i*</sub> is measured counter-clockwise from the north, *ρ*<sub>*i*</sub> is the angular separation of a component, and the centre of mass, (*x*<sub>*i*</sub>, *y*<sub>*i*</sub>) is the same in Cartesian coordinates. The instantaneous value of the argument of periastron is given as follows: ω(*t*) = ω₀ + ω(*t* − *T*<sub>p</sub>), where *T*<sub>p</sub> is the reference periastron epoch and ω₀ is the value of the periastron argument at the reference epoch. Instead of computing the semi-major axis for each component of a binary, the semi-major axis *a* and the mass ratio *q* = *M*<sub>1</sub>/*M*<sub>2</sub> were used; the semi-major axes of primary and secondary can be computed with the following formulae: *a*<sub>1</sub> = *aq*/(1 + *q*), *a*<sub>2</sub> = *a*/(1 + *q*). The periastron argument of the secondary is ω₂ = ω₁ + π.

In our application of Eqs. (5)–(8) component B is fixed at the beginning of the coordinate system because the observations are only sensitive to relative positions of the stars, not to the system as whole.

When the positions of all three components are known, objects representing each component can be placed at these positions. The uniform disk was chosen because all three components are detached and therefore only minor departures from spherical symmetry can be expected. The squared visibility V² and closure phase T₃φ for such a model can be computed analytically with the following formulae:

$$ |V_{S,k}(f)|^2 = \left| \frac{\sum_{j=1}^{N} L_{j,k} \frac{2J_1(\pi\theta_j B/\lambda_k)}{\pi\theta_j B/\lambda_k} e^{-2\pi i(f \cdot r)} |^2 \right. , \quad (9) $$

$$ T_3\phi_{S,k}(f_1, f_2) = \arg[V_{S,k}(f_1)V_{S,k}(f_2)V_{S,k}(-f_1 - f_2)], \quad (10) $$

where index *j* denotes a component of the triple system, *k* the spectral band, *V* the visibility, *f* = (*u*, *v*) the spatial frequency, *L* the luminosity fraction, *B* the length of the baseline, *θ* the diameter of the uniform disk, *λ* the effective wavelength (the central wavelength of the spectral band), *J*<sub>1</sub> the first-order Bessel function, *r* = (*x*<sub>*j*</sub>, *y*<sub>*j*</sub>) the Cartesian coordinates of a component computed with Eqs. (5)–(8), and *N* the total number of components in the system. The uniform disk diameter *θ* is also a wavelength-dependent quantity, therefore a different radius should be derived for each spectral band. Nonetheless, the dependency is very weak (order of 10⁻³ mas for the whole wavelength span of our data).

## 6.2. Orbital solution for all available spectro-interferometric observations

The model given by Eqs. (5)–(10) was fitted to calibrated squared visibilities from all four instruments, that is, CHARA/VEGA, NPOI, MARK III, and VLTI/AMBER. The best-fit set of parameters was determined using the least-squares

method, that is, by minimising the following χ²:

$$ \chi^2 = \sum_{k=1}^{N_F} \sum_{j=1}^{N_V} \left| \frac{|V_k|^2(f_j) - |V_{S,k}|^2(f_j)}{\sigma_k(f_j)} \right|^2 + \sum_{k=1}^{N_F} \sum_{j=1}^{N_T} \left| \frac{T_3\phi_k(f_{1,j}, f_{2,j}) - T_3\phi_{S,k}(f_{1,j}, f_{2,j})}{\sigma_k(f_{1,j}, f_{2,j})} \right|^2 \quad (11) $$

where V² (T₃φ) is the observed squared visibility (the observed closure phase), V<sub>S</sub>² (T₃φ<sub>S</sub>) the synthetic squared visibility computed with Eq. (9) (the synthetic closure phase, Eq. (10)), *f* = (*u*, *v*) the spatial frequency, *σ* the standard deviation of an observation, *N<sub>V</sub>* the total number of squared visibility observations, *N<sub>T</sub>* the total number of closure phase observations, and *N<sub>F</sub>* the total number of spectral bands.

The phase coverage of the inner and the outer orbits is good enough (see Fig. 1) to allow fitting of all orbital elements. Our strategy was to keep as many parameters free as possible, since this model is independent of those presented in Sects. 3 and 4. However, the angular size of the inner orbit is small and its ephemeris is obtained with greater precision by the photometry and spectroscopy. The eccentricity of orbit 1 was set to zero (see Table 5) because there were no signs of a significant eccentricity in previous analyses. A number of trial runs have shown that the inclination *i*<sub>1</sub> and the mass ratio *q*<sub>1</sub> are poorly constrained by the interferometric observations. If optimised, both converged to values not consistent with previous analyses (*i*<sub>1</sub> ≈ 78 ± 5 deg, *q*<sub>1</sub> = 0.8 ± 0.10). Investigation of χ² maps surrounding these values has shown large shallow valleys that spread up to regions with values consistent with photometric and spectroscopic models. To stay on the safe side, we fixed both parameters at values obtained from the spectroscopy and photometry because they were estimated with much higher precision.

The global minimum of Eq. (11) was determined with the differential evolution algorithm (Storn & Price 1997) and was locally optimised with the sequential least-squares routine (Kraft 1988). The parameters of the best-fitting model are listed in Table 11. A large portion of the parametric space was searched<fnref n="6" />. The initial parametric space was equally sampled with a population which consisted of 1500 members. The population evolved until the mean energy of the population (i.e. the mean χ² divided by its standard deviation and multiplied by the tolerance) was greater than one. The tolerance was set to 10⁻³ and the procedure took from 50 to 100 iterations to finish.

The final reduced χ²<sub>R</sub> ≈ 5.806 is much larger than 1 because the true uncertainty of the V² derived with the reduction pipeline is underestimated. The reason is that the high χ²<sub>R</sub> is given mainly by data that were acquired at low spatial resolution and are expected to be easiest to reduce. Another reason is that the angular slit width of all interferometric instruments is comparable to the angular separation of component C and the triple, meaning that it cannot be guaranteed that it was recorded. The full amplitude of squared visibility variations caused by component C ranges from ≤0.035 in the V band to ≤0.050 in the K band. It introduces systematic errors that we cannot correct for. The last reason are imperfections of the model. We had to accept several simplifications to stabilise the fit. Uncertainties of the best-fit

<footnote n="6">⁶ The investigated parametric space is given by the following ranges: θ<sub>B</sub> ∈ [0.0, 1.0] mas; L<sub>B</sub> ∈ [0.4, 0.8]; L<sub>Aa</sub> ∈ [0.1, 0.3]; T<sub>p,2</sub> ∈ [55 600.0, 55 620.0] RJD; a<sub>2</sub> ∈ [13, 18] mas; e<sub>2</sub> ∈ [0.1, 0.3]; i<sub>2</sub> ∈ [50, 130] deg; ω<sub>2</sub> ∈ [0, 180] deg; Ω<sub>2</sub> ∈ [0, 360] deg; ω̇<sub>2</sub> ∈ [1.5, 4.0] deg yr⁻¹; a<sub>1</sub> ∈ [1.0, 3.0] mas; Ω<sub>1</sub> ∈ [0, 360] deg.</footnote>
<table><caption>**Table 11.** Parameters corresponding to the best fit of all available interferometric observations with the model defined by Eqs. (5)–(10).</caption><thead><tr><th>Elements</th><th>Units</th><th colspan="3">Values</th></tr><tr><th>Component</th><th></th><th colspan="3">Component properties</th></tr><tr><th></th><th></th><th>B</th><th>Aa</th><th>Ab</th></tr></thead><tbody><tr><td>θ</td><td>(mas)</td><td>0.407 ± 0.031</td><td>0.252<sup>1.2</sup></td><td>0.231<sup>1.2</sup></td></tr><tr><td>L<sub>R</sub><sup>Δλ<sub>1</sub></sup></td><td></td><td>0.6373 ± 0.0085</td><td>0.197 ± 0.014</td><td>0.166 ± 0.016</td></tr><tr><td>L<sub>R</sub><sup>Δλ<sub>2</sub></sup></td><td></td><td>0.60<sup>1</sup></td><td>0.22<sup>1</sup></td><td>0.18<sup>1</sup></td></tr><tr><th>Orbit</th><th></th><th colspan="3">Orbital properties</th></tr><tr><th></th><th></th><th>2</th><th>1</th><th></th></tr><tr><td>P<sub>AN</sub></td><td>(d)</td><td>145.471 ± 0.045</td><td>—</td><td>—</td></tr><tr><td>P<sub>S</sub></td><td>(d)</td><td>145.150 ± 0.061</td><td>7.14664<sup>1</sup></td><td>—</td></tr><tr><td>T<sub>p</sub></td><td>(RJD)</td><td>55 609.36 ± 0.64</td><td>—</td><td>—</td></tr><tr><td>T<sub>min</sub></td><td></td><td>—</td><td>56 224.7248<sup>1</sup></td><td>—</td></tr><tr><td>a</td><td>(mas)</td><td>15.93 ± 0.070</td><td>1.89 ± 0.11</td><td>—</td></tr><tr><td>e</td><td></td><td>0.212 ± 0.0040</td><td>0.0<sup>1</sup></td><td>—</td></tr><tr><td>q</td><td></td><td>—</td><td>0.945<sup>1</sup></td><td>—</td></tr><tr><td>i</td><td>(deg)</td><td>86.67 ± 0.12</td><td>86.85<sup>1</sup></td><td>—</td></tr><tr><td>ω</td><td>(deg)</td><td>8.4 ± 1.6</td><td>90.0<sup>1</sup></td><td>—</td></tr><tr><td>Ω</td><td>(deg)</td><td>148.453 ± 0.066</td><td>148.435<sup>3</sup> ± 1.9</td><td>—</td></tr><tr><td>ω̇</td><td>(deg yr<sup>-1</sup>)</td><td>2.02 ± 0.31</td><td>0.0<sup>1</sup></td><td>—</td></tr><tr><td>χ<sub>R</sub><sup>2</sup></td><td></td><td></td><td>5.806</td><td>—</td></tr></tbody></table>

**Notes.** $\theta$ denotes the angular (uniform-disk) diameter. (1) The parameter was kept fixed. (2) Estimated from the solution from Table 9 and the HIPPARCOS parallax. (3) A solution shifted by 180 deg is also possible and has an identical reduced $\chi_R^2$. $\Delta\lambda_1$ = [500, 800] nm, and $\Delta\lambda_2$ = [1200, 2410] nm.

parameters were estimated at 68% confidence intervals from the $\chi_R^2$ scaled to one.

Several attempts have shown that we are insensitive to the diameters of components Aa and Ab, because we lack enough observations at very long baselines (reaching up to 300 m). If they were set free, the solution would converge to unrealistic values ($\ge$1.0 mas), therefore they had to be fixed at values given by the parallax of the system and the light-curve solution (see Table 9). Convergence of the orbital parameters of orbit 1 was in general slow because the bulk of observations (NPOI, AMBER) was taken at low spatial resolution, at which this orbit is almost mainly on observations from VEGA/CHARA.

Our model allows fitting separate sets of relative luminosities $L_R$ for each passband because the visibilities were estimated in narrow passbands: four for CHARA/VEGA, sixteen for NPOI, and ≈40 for VLTI/AMBER. It was not possible to divide the data into a larger number of small groups and to densely sample the relative luminosity of components Aa, Ab, and B as a function of the wavelength. After a set of trial attempts, we split the data into two subsets: visible (MARKIII, NPOI, CHARA/VEGA) and infrared (AMBER). This sampling is justified by the very low variability of the luminosity ratios with wavelength of all stars within the visible and infrared regions, which we checked using synthetic spectra from the PHOENIX grid (Husser et al. 2013). The relative luminosities of components Aa and Ab did not converge to plausible values for the infrared subset (it generally predicted a too low luminosity ratio between the two components of orbit 1), therefore we decided to use the estimate based on the PHOENIX grid and radii obtained from the light-curve analysis for components Aa and Ab, and the radius of component B was taken from Harmanec (1988).

The best-fit set of parameters is listed in Table 11 and a plot of the model vs. the observations is shown in Figs. C.1–C.10. The model qualitatively fits the variations of the $V^2$ (i.e. the

curvature of the model data agrees with the curvature of the observed $V^2$) for all spectro-interferometric data very well.

# 7. Summary of analyses based on simple analytic models

Here we critically compare the results of individual observational methods and derive the properties of the system.

## 7.1. Performance of different observational methods

Despite the subtitle, the individual models we used to evaluate different observational methods were not completely independent because the results from one method often served as a starting point for another. In some cases it was mandatory to take a parameter value from another model to stabilise the convergence to a steady solution. In the following paragraphs we discuss the outcome of different methods and their accuracy. An overview of all fitted parameters is given in Table 12 obtained through different methods (i.e. more values are given for some parameters). Corresponding properties of the orbits and stars are also listed. Orbital elements of orbit 3 are not listed because their properties were constrained only by astrometry, and they are presented separately in Table 10. The mass of component C is briefly discussed here.

- *The spectroscopic elements:* elements (K, e, T<sub>p</sub>, P, ω, ω̇) of both orbits are estimated better from the fit of directly measured RVs with an analytic model (see Table 5, Eqs. (2) and (3)). The spectral disentangling works with a much more complex model, and the resulting orbital elements depend on the shape of the separated profiles (and vice versa), which come out warped (the degree of the warp is shown by grey line in Fig. A.1). The warp is most pronounced for
<table><caption>**Table 12.** Summary of parameters derived from the spectroscopic, photometric, and spectro-interferometric analyses.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Source</th><th colspan="2">Value</th></tr><tr><td></td><td></td><td></td><td colspan="2" style="text-align: center;">Orbital properties</td></tr><tr><td></td><td></td><td></td><td>2</td><td>1</td></tr></thead><tbody><tr><td>P<sub>AN</sub></td><td>(d)</td><td>RV/LC</td><td>145.579 ± 0.048</td><td>7.14664 ± 0.00002</td></tr><tr><td>P<sub>S</sub></td><td>(d)</td><td>RV/LC</td><td>145.113 ± 0.071</td><td>7.14664 ± 0.00002</td></tr><tr><td>T<sub>p</sub></td><td>(RJD)</td><td>RV/-</td><td>55 609.46 ± 0.52</td><td>-</td></tr><tr><td>T<sub>min</sub></td><td>(RJD)</td><td>-/LC</td><td>-</td><td>56 224.72482 ± 0.00022</td></tr><tr><td>a</td><td>(R<sub>☉</sub>)</td><td>IF+HP/IF+HP</td><td>219 ± 15</td><td>26.1 ± 2.3</td></tr><tr><td>a</td><td>(R<sub>☉</sub>)</td><td>IF+RV/IF+RV</td><td>* 229.0 ± 7.7</td><td>-</td></tr><tr><td>a</td><td>(R<sub>☉</sub>)</td><td>-/RV+LC</td><td>-</td><td>* 25.550 ± 0.090</td></tr><tr><td>a<sub>angular</sub></td><td>(mas)</td><td>IF/IF</td><td>15.93 ± 0.10</td><td>1.89 ± 0.11</td></tr><tr><td>e</td><td></td><td>RV/-</td><td>* 0.2101 ± 0.0053</td><td>0.0<sup>1</sup></td></tr><tr><td>e</td><td></td><td>IF/-</td><td>0.2120 ± 0.0040</td><td>0.0<sup>1</sup></td></tr><tr><td>q</td><td></td><td>RV/RV</td><td>0.889 ± 0.056</td><td>0.9438 ± 0.0036</td></tr><tr><td>i</td><td>(deg)</td><td>IF/IF</td><td>86.67 ± 0.12</td><td>86.85<sup>1</sup></td></tr><tr><td>i</td><td>(deg)</td><td>-/LC</td><td>-</td><td>* 86.85 ± 0.22</td></tr><tr><td>ω</td><td>(deg)</td><td>RV/-</td><td>* 9.25 ± 1.42</td><td>90.0<sup>1</sup></td></tr><tr><td>ω</td><td>(deg)</td><td>IF/-</td><td>8.4 ± 1.6</td><td>90.0<sup>1</sup></td></tr><tr><td>ω</td><td>(deg yr<sup>-1</sup>)</td><td>RV/-</td><td>2.90 ± 0.33</td><td>0.0<sup>1</sup></td></tr><tr><td>ω</td><td>(deg yr<sup>-1</sup>)</td><td>IF/-</td><td>2.02 ± 0.31</td><td>0.0<sup>1</sup></td></tr><tr><td>Ω</td><td>(deg)</td><td>IF/IF</td><td>148.453 ± 0.066</td><td>148.4<sup>2</sup> ± 1.9</td></tr><tr><td></td><td></td><td></td><td colspan="2" style="text-align: center;">Component properties</td></tr><tr><td>Component</td><td></td><td></td><td>B</td><td>Aa</td><td>Ab</td></tr><tr><td>T<sub>eff</sub></td><td>(K)</td><td>SP</td><td>14 190 ± 150</td><td>10 700 ± 160</td><td>* 10 480 ± 130</td></tr><tr><td>log g<sub>[cgs]</sub></td><td>(K)</td><td>LC</td><td>-</td><td>10700<sup>1</sup></td><td>10450 ± 150</td></tr><tr><td></td><td></td><td>SP</td><td>4.527 ± 0.041</td><td>4.08 ± 0.12</td><td>4.01 ± 0.10</td></tr><tr><td></td><td></td><td>RV+LC</td><td>4.09 ± 0.11</td><td>4.330 ± 0.019</td><td>4.348 ± 0.022</td></tr><tr><td>ν sin i</td><td>(km s<sup>-1</sup>)</td><td>SP</td><td>229.2 ± 1.7</td><td>12.6 ± 2.6</td><td>14.3 ± 3.1</td></tr><tr><td>m</td><td>(M<sub>☉</sub>)</td><td>RV+LC</td><td>* 3.89 ± 0.25</td><td>* 2.252 ± 0.027</td><td>* 2.125 ± 0.027</td></tr><tr><td>m</td><td>(M<sub>☉</sub>)</td><td>RV+IF</td><td>3.60 ± 0.52</td><td>2.08 ± 0.48</td><td>1.96 ± 0.45</td></tr><tr><td>R</td><td>(R<sub>☉</sub>)</td><td>RV+LC</td><td>-</td><td>1.700 ± 0.035</td><td>1.618 ± 0.039</td></tr><tr><td>R</td><td>(R<sub>☉</sub>)</td><td>IF+HP</td><td>2.81 ± 0.28</td><td>-</td><td>-</td></tr><tr><td>θ</td><td>(mas)</td><td>IF</td><td>0.407 ± 0.031</td><td>-</td><td>-</td></tr><tr><td>M<sub>BOL</sub></td><td>(mag)</td><td>LC+HP</td><td>-</td><td>0.247 ± 0.017</td><td>0.235 ± 0.017</td></tr><tr><td>V</td><td>(mag)</td><td>LC+RV+IF</td><td>-1.14 ± 0.22</td><td>0.923 ± 0.079</td><td>1.120 ± 0.075</td></tr><tr><td>B - V</td><td>(mag)</td><td>LC</td><td>4.250 ± 0.10</td><td>5.46 ± 0.11</td><td>5.63 ± 0.11</td></tr><tr><td>U - B</td><td>(mag)</td><td>LC</td><td>-0.12 ± 0.16</td><td>-0.05 ± 0.16</td><td>-0.03 ± 0.14</td></tr><tr><td>U - B</td><td>(mag)</td><td>LC</td><td>-0.446 ± 0.16</td><td>-0.09 ± 0.14</td><td>-0.07 ± 0.14</td></tr><tr><td>V<sub>0</sub></td><td>(mag)</td><td>LC</td><td>4.24 ± 0.65</td><td>5.54 ± 0.65</td><td>5.68 ± 0.65</td></tr><tr><td>B<sub>0</sub> - V<sub>0</sub></td><td>(mag)</td><td>LC</td><td>-0.120 ± 0.085</td><td>-0.018 ± 0.167</td><td>-0.015 ± 0.162</td></tr><tr><td></td><td></td><td></td><td colspan="3" style="text-align: center;">Parallax</td></tr><tr><td>π<sub>a1</sub></td><td>(mas)</td><td></td><td colspan="3">15.91 ± 0.93</td></tr><tr><td>π<sub>a2</sub></td><td>(mas)</td><td></td><td colspan="3">* 14.96 ± 0.51</td></tr><tr><td>π<sub>DM,Aa</sub></td><td>(mas)</td><td></td><td colspan="3">14.3 ± 4.3</td></tr><tr><td>π<sub>DM,Ab</sub></td><td>(mas)</td><td></td><td colspan="3">14.4 ± 4.4</td></tr><tr><td>π<sub>DM,B</sub></td><td>(mas)</td><td></td><td colspan="3">13.3 ± 2.2</td></tr></tbody></table>

**Notes.** In some cases more values are listed for a parameter, to show that the methods do not contradict each other. Because they were safely resolved only with the astrometry, elements of orbit 3 are not listed here, but in Table 10 and the mass of component C is estimated and briefly discussed in Sect. 7. The listed parameters are the anomalous period $P_{AN}$, the sidereal period $P_{S}$, the periastron epoch $T_{p}$, the epoch of the primary minimum $T_{min}$, the semi-major axis $a$, the mass ratio $q$, the eccentricity $e$, the inclination $i$, the periastron argument $ω$, the position angle of the nodal line $Ω$, the effective temperature $T_{eff}$, the surface gravitational acceleration $log g$, the projected rotational velocity $ν sin i$, the mass $m$, the radius $R$, the angular diameter $θ$, the bolometric magnitude $M_{BOL}$, the Johnson V magnitude and colour indices $V, B - V, U - B$, the dereddened Johnson V magnitude and colour index $V_{0}, B_{0} - V_{0}$, and the parallax $π$. (1) Assumed. (2) A solution where $Ω_{1} = 328.4 ± 1.9$, is also plausible and has identical $χ^{2}$. (3) Parameters that are likely the closest to the true nature of $ξ$ Tau. Sources: RV solution of the RV curve presented in Table 5, SP comparison of the observed and synthetic spectra presented in Table 7, LC solution of the light curve presented in Table 9, IF solution of the $V^{2}$ and $T_{3}φ$ presented in Table 11, HP. the HIPPARCOS parallax $π = 15.60 ± 1.04$ mas. The parallaxes: $π_{a1}$ estimated from the size of the semi-major axis of orbit 1 (physical and angular), $π_{a2}$ estimated from the size of the semi-major axis of orbit 2 (physical and angular), $π_{DM,Aa}$ estimated from the distance modulus of component Aa, $π_{DM,Ab}$ estimated from the distance modulus of component Ab, $π_{DM,Ab}$ estimated from the distance modulus of component B.

component B, meaning that especially the mass ratio $q_{2}$ coming from the method cannot be trusted. On the other hand, the thin lines of components Aa and Ab constrain the RVs very well even if the separated spectrum is not perfect, and for the

remaining orbital parameters the disentangling therefore provides values that agree with the fit of directly measured RVs. – *The ephemeris of orbit 1:* the photometric solution presented in Table 9 yields the best ephemeris ($T_{min,1}, P_{1}$) of orbit 1
especially thanks to high-precision observations from the satellite MOST. The ephemeris for orbit 1 estimated from the RVs does not agree within uncertainties with the photometric one. It can be caused by the lower precision of RV measurements around eclipses.

- *The eccentricity of orbit 1:* it was set to zero throughout the analyses because the precision of data does not allow a reliable determination. The analysis of the light curve from the satellite MOST shows a hint of a small eccentricity, but the relative position of minima is also affected by ETVs, and we are unable to discern one from the other with the analytic models. The dynamics of the system (see Sects. 8 and 9) shows that the eccentricity should oscillate with an amplitude $\Delta e \approx 0.01$. This introduces a jitter of the relative position of the primary and secondary minimum and increases uncertainty of the radii when a circular model is applied.

- *The inclination of orbit 1:* it is determined accurately using the light-curve analysis presented in Table 9. The value obtained from the interferometric model suffers from large uncertainty and is about 10 deg off the photometric solution. This is probably caused by the low number of observations at high spatial frequencies and the calibration systematic errors, which are likely more pronounced for high-frequency data.

- *The longitude of the ascending node:* the longitude of the ascending node of orbit 1 has a mirror solution $\Omega_1 = \Omega_1 + 180$ deg with (almost) the same value of the $\chi^2_R$, while the $\Omega_2$ is determined uniquely because the NPOI instrument acquired a large number of closure phase measurements. This means that it is not possible to say whether the motion of orbit 1 relative to orbit 2 is prograde or retrograde based solely on the spectro-interferometric data.

- *The relative luminosities:* they were determined from the light-curve solution, the comparison of synthetic and observed spectra, and from the interferometric solution.

* The light-curve solution best describes their variations with the wavelength, but the values suffer from large uncertainties because of correlations between the fitted parameters.

* The fit of synthetic spectra to observed ones is quite insensitive to relative luminosities, but this is the case only because small parts of red spectra were fitted that contain only three weak spectral lines. The relative luminosities obtained in the regions around $H\gamma$ and $H\beta$ roughly agree with the values obtained for the B band from the light-curve solution.

* The bulk of the interferometric observations falls somewhere between the V and R bands. Therefore the relative luminosities detected with the spectro-interferometry are close to the V-band value obtained from the light-curve solution. We were not able to obtain plausible estimates of relative luminosities for the infrared subset (AMBER) because the observations have low spatial resolution and do not resolve the eclipsing binary well.

- *The effective temperatures:* they are given better by the fits of observed spectra to synthetic ones because the fitted regions contain many spectral lines (especially the region $\Delta\lambda = [4280, 4495]$ Å) where the photometry relies on four broad-band filters alone. In addition, Prša & Zwitter (2006) stated that it is not possible to obtain accurate effective temperatures of the two components of an eclipsing binary from the light-curve solution unless the colour-constraining method (described by them) is employed. According to the authors, the problem is even more pronounced when the two

components are alike. Therefore we fixed the primary temperature and only optimised the secondary temperature. The result agrees with that obtained from the comparison of observed and synthetic profiles within the respective errors. The spectral types corresponding to these temperatures are B9 for components Aa and Ab and B5-6 for component B.

- *The semi-major axes and masses:* the physical size of the semi-major axes derived from the spectro-interferometry and the HIPPARCOS parallax (orbits 1 and 2) and those derived from the spectroscopy and photometry (orbit 1) and spectroscopy and spectro-interferometry (orbit 2) agree with each other within their uncertainties. The same applies to masses, which seem to fall within the limits of normal main-sequence (MS hereafter) masses corresponding to the respective spectral types (Harmanec 1988) – $m^{Aa} = 2.25 \pm 0.03 \in [1.71, 2.41]$ $M_\odot$, $m^{Ab} = 2.13 \pm 0.03 \in [1.71, 2.41]$ $M_\odot$, $m^B = 3.89 \pm 0.25 \in [3.63, 4.6]$ $M_\odot$.

- *The total mass of the system and mass of component C:* using the parallax $\pi_{a2} = 14.96 \pm 0.51$ and the solution presented in Table 10, we can estimate the total mass of the system $m^{Aa+Ab+B+C} = 9.88 \pm 1.06$ $M_\odot$. A comparison with the masses of the inner triple subsystem gives an estimate of the mass of component C $m^C = 1.61 \pm 1.18$ $M_\odot$ that agrees with early F-type or late A-type star.

- *The component radii:* all components seem to have normal radii for their respective spectral type (again checked against Harmanec 1988) – $R^{Aa} = 1.70 \pm 0.04 \in [1.40, 2.06]$ $R_\odot$, $R^{Ab} = 1.62 \pm 0.04 \in [1.40, 2.06]$ $R_\odot$, $R^B = 2.8 \pm 0.3 \in [2.13, 2.85]$ $R_\odot$.

- *The dereddened colour index B-V:* these are derived with a high level of uncertainty because of the high uncertainty in the luminosity ratios in different bands and the uncertainty of bolometric magnitudes. We compared the dereddened colour indices against tables computed by Flower (1996),

$$T_{\text{teff,FLOWER}}^{\text{B}}(-0.120) = 12\,370 \text{ K},$$

$$T_{\text{teff,FLOWER}}^{\text{Aa}}(-0.018) = 9810 \text{ K}, \text{ and}$$

$$T_{\text{teff,FLOWER}}^{\text{Aa}}(-0.015) = 9760 \text{ K}.$$

They very roughly agree with the values found by the comparison of the observed and synthetic spectra. The uncertainty bars of the colour indices are very generous and match a wide range of temperatures.

- *The distance:* the number of applied observational methods allows us to estimate the distance of $\xi$ Tau from the ratio of the physical and angular size of the semimajor axes and from the distance modulus. The former seems to prefer parallax, which is slightly lower than the HIPPARCOS parallax (but still within error bars), the latter also places $\xi$ Tau farther than the HIPPARCOS observations, but their uncertainties are large, meaning that they do not contradict the HIPPARCOS parallax. The parallax estimated from the ratio of the physical and angular size of the semi-major axis of the outer orbit yields the most precise parallax, $\pi_{a2} = 14.96 \pm 0.51$ mas.

## 7.2. Conclusion of the analytic models

The spectroscopy, the photometry, and the interferometry were studied with traditional (semi-) analytic models. We found that results obtained from different methods are consistent with each other, although some of them give better estimates of a particular set of parameters than others. We took advantage of this differential sensitivity and compiled a resulting set of fundamental properties of the system.
During the analyses described in previous sections, we noted two effects that indicate the dynamical interaction in ξ Tau: the advance of the apsidal line of orbit 2, and the eclipse timing variations (ETVs) in system 1. The first effect was explicitly taken into account because omitting it would cause significant inconsistency between observations and model. The latter effect was almost overlooked if it had not been for the indication in the very accurate photometric data from the MOST satellite. However, the analytic models above give only limited insights into dynamical effects in a four-body system such as ξ Tau. Nonetheless, they provide very good results that are also needed as a starting point for a more sophisticated solution based on an approach that includes dynamical evolution in a more complete way. We proceed in two steps.

In Sect. 8 we develop a numerical model that consistently takes into account the gravitational interaction of all stars in the ξ Tau. We use a fully numerical implementation, basically a standard N-body integrator, which we extend by subroutines that allow us to model several types of observables relevant for the ξ Tau dataset.

Next, in Sect. 9 we summarise relevant analytic formulae obtained by methods of perturbation theory, which provide insights into results from the fully numerical approach in Sect. 8. Despite their limitations, we find the analytic formulation of the most important orbital perturbations useful. It does not only allow us to understand basic features in the numerical integrations, but also readily provides the parametric dependencies.

## 8. N-body model of ξ Tauri with mutual interactions

The quadruple nature of ξ Tauri and its relatively compact packing require us to proceed with an advanced N-body model that can account for mutual gravitational interactions of all four components. To this point, we now describe our numerical integrator, a definition of a suitable $\chi^2$ metric, and the overall results of our fitting procedure.

### 8.1. Numerical integrator and $\chi^2$ metric

We use a standard Bulirsch–Stoer N-body numerical integrator from the SWIFT package (Levison & Duncan 2013). Our method is quite general. We can model classical Keplerian orbits, of course, but also non-Keplerian orbits (involving N-body interactions). We treat all stars as point masses only, however. We have no higher-order gravitational terms and no tides in our model.

As explained below, this is a significant improvement of our previous application in Brož et al. (2010) because we can now account not only for the light-time effect, but for complete eclipse timing variations (ETVs) of the inner binary that arise from both direct and indirect gravitational perturbations. At the same time, we do not use the simplification of Brož et al. (2010) and consider all the components separately because the equivalent gravitational moment

$$J_2 \approx \frac{1}{2} \left(\frac{a_1}{r}\right)^2 \frac{m^{\text{Aa}} m^{\text{Ab}}}{(m^{\text{Aa}} + m^{\text{Ab}})^2} \approx 2 \times 10^{-3} \quad (12)$$

of the inner eclipsing binary Aa+Ab is large at the distance of the component B.

We used five different coordinate systems: (i) Aa-centric (to generally specify initial conditions and eclipse detection); (ii) barycentric (for the numerical integration itself); (iii) Aa+Ab photocentric (to compare with interferometric observations of

component B); (iv) Aa+Ab+B photocentric (ditto for component C); and (v) Jacobian (to compute hierarchical orbital elements).

Initial conditions at a given epoch $T_0$ can be specified either in Cartesian coordinates with $x, y$ in the sky plane and $z$ in the radial direction, or in osculating orbital elements. This very choice has a substantial role because the outcome of the fitting procedure will be generally (slightly) different. The orbital elements can be considered less strongly correlated quantities than Aa-centric Cartesian coordinates.

We accounted for as many observational data as possible using the following joint metric<fnref n="7" />:

$$\chi^2 = \chi_{\text{rv}}^2 + \chi_{\text{etv}}^2 + \chi_{\text{edv}}^2 + \chi_{\text{sky}}^2, \quad (13)$$

$$\chi_{\text{rv}}^2 = \sum_{j=1}^{4} \sum_{i=1}^{N_{\text{rv},j}} \frac{(v'_{\text{zb},ji} + \gamma - v_{\text{rad},ji})^2}{\sigma_{\text{rv},ji}^2}, \quad (14)$$

$$\chi_{\text{etv}}^2 = \sum_{i=1}^{N_{\text{etv}}} \frac{(t'_{Ai} - t_{Ai})^2}{\sigma_{\text{etv},i}^2}, \quad (15)$$

$$\chi_{\text{edv}}^2 = \sum_{i=1}^{N_{\text{edv}}} \frac{(\epsilon'_{Ai} - \epsilon_{Ai})^2}{\sigma_{\text{edv},i}^2}, \quad (16)$$

$$(\Delta x_{ji}, \Delta y_{ji}) = \mathbf{R} \left(-\phi_{\text{ellipse}} - \frac{\pi}{2}\right) \times \begin{pmatrix} x'_{p,ji} - x_{p,ji} \\ y'_{p,ji} - y_{p,ji} \end{pmatrix}, \quad (17)$$

$$\chi_{\text{sky}}^2 = \sum_{j=3}^{4} \sum_{i=1}^{N_{\text{sky},j}} \left\{ \frac{(\Delta x_{ji})^2}{\sigma_{\text{sky major},ji}^2} + \frac{(\Delta y_{ji})^2}{\sigma_{\text{sky minor},ji}^2} \right\}, \quad (18)$$

where the notation is briefly described in Table 13. The dashed quantities are the model values linearly interpolated to the exact times $t_i$ of observations. The index $j$ goes over the list of components Aa, Ab, B, C (i.e. $j = 1 = Aa, \dots$), while the index $i$ corresponds to the observational data.

In our N-body model we do not fit the observed spectra using synthetic ones, individual light curve points, or interferometric fringes. We use higher-level observational data instead that were reduced and derived in previous sections. Hence we fit RV measurements for the three components Aa, Ab, and B, altogether $N_{\text{rv}} = 843$, minima timings for the eclipses in the inner binary (Aa+Ab), $N_{\text{etv}} = 35$, and astrometric observations for components B and C, $N_{\text{sky}} = 49$. The latter is a subset of measurements from NPOI and WDS, for which it was possible to convert fringe visibilities (averaged over one night) into distance–angle values. The individual uncertainties of the observations used in this section were modified as follows: $\sigma_{\text{rv}} \ge 2 \text{ km s}^{-1}$ due to calibration uncertainties, $\sigma_{\text{etv}} \ge 0.001 \text{ d} = 1.5 \text{ min}$ because the quasi-periodic oscillations visible in the MOST light curve shift minima timings in a random fashion, and $\sigma_{\text{sky}} = 3 \text{ mas}$ (as in Tokovinin et al. 2015) or 5 mas if not reported in WDS.

We assumed the nominal distance $d = 64.1 \text{ pc}$ for ξ Tau. The stellar radii for an eclipse detection were $R^{\text{Aa}} = 1.700 R_{\odot}$ and

<footnote n="7">⁷ The program used for these computations, including sources and all input data, is available at http://sirrah.troja.mff.cuni.cz/~mira/xitau/</footnote>
<table><caption>**Table 13.** Notation used for various coordinates, velocities, and uncertainties that we used in our N-body model.</caption><tr><td>x, y, z</td><td>Aa-centric coordinates</td></tr><tr><td>v<sub>x</sub>, v<sub>y</sub>, v<sub>z</sub></td><td>Aa-centric velocities</td></tr><tr><td>x<sub>pB</sub>, y<sub>pB</sub></td><td>Aa+Ab photocentric sky-plane coordinates<br>(of component B)</td></tr><tr><td>x<sub>pC</sub>, y<sub>pC</sub></td><td>Aa+Ab+B photocentric coordinates<br>(of component C)</td></tr><tr><td>v<sub>zb</sub></td><td>barycentric radial velocity</td></tr><tr><td>γ</td><td>systemic velocity</td></tr><tr><td>v<sub>rad</sub></td><td>observed radial velocity</td></tr><tr><td>t<sub>A</sub></td><td>mid-epoch of an eclipse of the Aa+Ab pair</td></tr><tr><td>ε<sub>A</sub></td><td>eclipse duration</td></tr><tr><td>σ<sub>tv</sub></td><td>uncertainty of the radial velocity</td></tr><tr><td>σ<sub>etv</sub></td><td>uncertainty of the eclipse mid-epoch timing</td></tr><tr><td>σ<sub>edv</sub></td><td>uncertainty of the eclipse duration</td></tr><tr><td>σ<sub>sky major, minor</sub></td><td>uncertainty of the astrometric position,<br>angular sizes of the uncertainty ellipse</td></tr><tr><td>φ<sub>ellipse</sub></td><td>position angle of the ellipse</td></tr><tr><td>R(...</td><td>the corresponding 2 × 2 rotation matrix</td></tr></table>

$R^{\text{Ab}} = 1.612 R_{\odot}$, in agreement with the photometric inversion. The expected correlation among $R^{\text{Aa}}$, $R^{\text{Ab}}$, eclipse depth, eclipse duration and third light contribution is removed to some extent through spectroscopic observations (cf. Table 9).

The synthetic minimum distance $\Delta'$ between components Aa and Ab in the sky plane was determined analytically as the distance of the piece-wise straight line ($x_{\text{Ab}}$, $y_{\text{Ab}}$) from the origin in the Aa-centric coordinates, as provided by the numerical integration. The condition for an eclipse is then $\Delta' \le R^{\text{Aa}} + R^{\text{Ab}}$ and the corresponding time $t'_A$ is linearly interpolated from neighbouring points. The eclipse duration is then given by a simple geometry, $\epsilon'_A = 2 \sqrt{(R^{\text{Aa}} + R^{\text{Ab}})^2 - \Delta'^2/\bar{v}$, where $\bar{v}$ denotes the average velocity between the points. We thus straightforwardly account for disappearing eclipses and their durations, but we do not model (possible) eclipse depth variations at this stage.

To remove minor systematics in minima timings and eclipse duration, we attempted to suppress quasi-periodic oscillations visible in the MOST light curve by subtracting a function of the following form:

$$f(t) = C_0 + C_1(t - T_1) - [A_0 + A_1(t - T_1)] \sin \left[ \frac{2\pi(t - T_1)}{P_0 + P_1(t - T_1)} \right]. \quad (19)$$

Its coefficients ($C_0, C_1, T_1, A_0, A_1, P_0, P_1$) were always determined by a local fit in the surroundings of the given minimum. The resulting data are reported in Table 14.

The relative luminosities for photocentre computations were set to $L^{\text{Aa}} = 0.204$, $L^{\text{Ab}} = 0.174$, and $L^{\text{B}} = 0.622$, again in agreement with photometric observations.

Mass constraints also arise from the spectroscopic classification of the $\xi$ Tau components (A9 V, A9 V, B5 V, and F V). We can easily enforce reasonable limits for the component masses with the following artificial term:

$$\chi^2_{\text{mass}} = \sum_{j=1}^{4} \left[ \left( m_j - \frac{m_{j \text{min}} + m_{j \text{max}}}{2} \right) \frac{2}{m_{j \text{max}} - m_{j \text{min}}} \right]^{100}, \quad (20)$$

where we used $m^{\text{Aa}}$ and $m^{\text{Ab}} \in (0.9, 3.0) M_{\odot}$, $m^{\text{B}} \in (3.5, 3.9) M_{\odot}$, $m^{\text{C}} \in (0.9, 2.0) M_{\odot}$ as the limits; the exponent is rather arbitrary.

The integrator and its internal time step were controlled by the parameter $\epsilon_{\text{BS}} = 10^{-8}$ (unitless), which ensures a sufficient

<table><caption>**Table 14.** Subset of minima timings $t_A$ and eclipse durations $\epsilon_A$ determined from MOST light curves, which were corrected for quasi-periodic oscillations by means of Eq. (19), and corresponding uncertainties $\sigma_{\text{etv}}$ and $\sigma_{\text{edv}}$.</caption><thead><tr><th>t<sub>A</sub><br>RJD</th><th>σ<sub>etv</sub><br>day</th><th>ε<sub>A</sub><br>day</th><th>σ<sub>edv</sub><br>day</th></tr></thead><tbody><tr><td>56 224.7242</td><td>0.0010</td><td>0.2656</td><td>0.0069</td></tr><tr><td>56 228.3017</td><td>0.0012</td><td>0.2611</td><td>0.0035</td></tr><tr><td>56 231.8686</td><td>0.0010</td><td>0.2678</td><td>0.0069</td></tr><tr><td>56 235.4452</td><td>0.0010</td><td>0.2573</td><td>0.0035</td></tr></tbody></table>

accuracy. The integration time span was 1000 d forward and 11 000 d backward, and the output timestep $\Delta t = 0.5$ d for initial runs. We verified that this sampling is sufficient even for the trajectory with the strongest curvature and all necessary interpolations to the times of observations. For the final optimisation we decreased the value further to $\Delta t = 0.1$ d to suppress interpolation errors.

We used a standard simplex algorithm (Press et al. 1993) to search for local minima of $\chi^2$. We have 23 potentially free parameters, masses $m_j$, coordinates $x_j, y_j, z_j$, velocities $v_{xj}, v_{yj}, v_{zj}$ in the Aa-centric frame, or, alternatively, masses $m_j$ and three sets of orbital elements $a_j, e_j, I_j, \Omega_j, \omega_j, M_j$ in Jacobian coordinates, and the systemic velocity $\gamma$. The convergence tolerance for $\chi^2$ was set to $\epsilon_{\text{tol}} = 10^{-6}$, and the maximum number of iterations to 10 000 or to as low as 300 for extended surveys of the parameter space. We verified that this low number is sufficient to quickly detect local minima or to exclude their existence.

The initial epoch $T_0 = 2456 224.724705$ is very close to the first precise minimum of the MOST light curve. We can thus (almost) fix $x_{\text{Ab}} \approx y_{\text{Ab}} = 0$. At the same time, it is possible to (approximately) fix positions $x_{\text{pB}}, y_{\text{pB}}$ and $x_{\text{pC}}, y_{\text{pC}}$, derived by interferometry for an epoch close to $T_0$.

## 8.2. Resulting best fits

As expected, the 23-dimensional parameter space is vast and full of local minima, even at high $\chi^2$. We proceeded sequentially to avoid complications and used a set with 2012 data only, a set with data from 2011–2013, and one set with all observational data. Next we performed a survey of the parameter space (to ensure we did not miss an obvious global minimum), an optimisation of individual orbits (2 and 3), the mutual inclination of orbits 1 and 2, and then we switched from Cartesian coordinates to orbital elements. Finally, we let all parameters be free. The optimisation means that we started the simplex from scratch many times (with different initialisation) and let it converge (for a limited number of iterations). Our largest survey consisted of $10^5$ simplex runs, 300 steps each, that is, $3 \times 10^7$ models in total, so that we are confident that there is no other hidden minimum, at least within the ranges searched so far$^8$.

We are aware of three mirror solutions (and $2^3$ combinations), namely the inner binary can orbit in a retrograde or prograde sense with respect to orbit 2, so that $i'_1 = 180^\circ - i_1$. Moreover, its node can be shifted by $180^\circ$, $\Omega'_1 = \Omega_1 + 180^\circ$. Last

<footnote>$^8$ The ranges expressed in Cartesian coordinates were $z_{\text{Ab}} \in (-0.148, -0.088)$ au, $z_{\text{B}} \in (-1.47, -0.87)$ au, $z_{\text{C}} \in (-8.72, -2.72)$ au, $v_{x\text{Ab}} \in (-0.092, -0.032)$ au d$^{-1}$, $v_{y\text{Ab}} \in (0.050, 0.110)$ au d$^{-1}$, $v_{x\text{B}} \in (-0.078, -0.018)$ au d$^{-1}$, $v_{y\text{B}} \in (0.042, 0.102)$ au d$^{-1}$, $v_{z\text{B}} \in (-0.022, 0.038)$ au d$^{-1}$, $v_{x\text{C}} \in (-0.082, -0.022)$ au d$^{-1}$, $v_{y\text{C}} \in (0.025, 0.085)$ au d$^{-1}$, and $v_{z\text{C}} \in (-0.030, 0.030)$ au d$^{-1}$.</footnote>
<figure><img src="image_10.png" /></figure>

**Fig. 9.** One of the best-fit solutions for the ξ Tau system with our N-body model and using all available observational data. In this case, the resulting total χ² is 2578, while the number of degrees of freedom *v* = 908. **Top:** radial velocities *v*<sub>zbAa</sub>, *v*<sub>zbAb</sub>, *v*<sub>zbB</sub>, *v*<sub>zbc</sub> of the individual components; model values are denoted by lines (component Aa is black, not clearly visible, Ab grey, B blue, and C orange), observations by black error bars and residuals by thick red lines. **Middle:** O–C values for both primary and secondary minima timings; model timings are denoted by black points (very densely packed), observations by grey crosses, and O–C with its uncertainty by red error bars. **Bottom left:** astrometric positions of component B based on NPOI interferometric observations; model orbit *x*<sub>pB</sub>, *y*<sub>pB</sub> with respect to photocentre Aa+Ab (i.e. not w.r.t. B, as usually) is again denoted by a blue line, observations by black error bars and residuals by thick red lines. The orbit is not a single ellipse, but rather a complex trajectory that quickly precesses and is moreover affected by (slight) photocentre motions. **Bottom right:** similarly, astrometric positions of the distant component C *x*<sub>pC</sub>, *y*<sub>pC</sub> with respect to the Aa+Ab+B photocentre is denoted by an orange line. Component B is relatively luminous, which makes the orbit in these photocentric coordinates slightly jagged.

but not least, orbit 3 can have the opposite inclination, $i'_3 = -i_3$
(we have no direct RV measurements). These ambiguities are
discussed and partly resolved in the following paragraphs.

Our best fit is presented in Fig. 9 and Table 15. We note
that this is not the only fit that seems reasonable; there are many

more available in the surroundings. This can be partially seen in
Fig. 10 where one-dimensional χ² maps exhibit relatively broad
minima for the plotted parameters. Consequently, if we were to
use simplex within these ranges, we would surely find a different
minimum with slightly larger χ² (or even slightly smaller).
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Example of the one-dimensional $\chi^2$ mapping used to derive uncertainties of orbital elements for the $\xi$ Tau system. The dependencies of the $\chi^2$ values on the three nodes $\Omega_1, \Omega_2, \Omega_3$ are shown, while the remaining elements correspond to the best-fit values from Table 15. The preferred solution for $\Omega_1 \approx 331^\circ$, with $\chi^2 = 2578$, and a hint of a mirror solution at $\Omega_1 \approx 151^\circ$ are clearly visible. If the latter is optimised separately, we would obtain $\chi^2$ as low as 2749. The sudden increase of $\chi^2_{etv}$ corresponds to the disappearance of the eclipses of the inner binary, which naturally results in extreme O–C’s.</figcaption></figure>

<table><caption>**Table 15.** Initial osculating orbital elements $a_j, e_j, i_j, \Omega_j, \omega_j, M_j$ of the $\xi$ Tau system as derived by our N-body model.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th><th>Note</th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>0.1175673<sup>+0.0000007</sup><sub>-0.0000007</sub></td><td>a<sub>2</sub></td><td>1.08296<sup>+0.00033</sup><sub>-0.00031</sub></td><td>a<sub>3</sub></td><td>28.35<sup>+0.81</sup><sub>-0.78</sub></td><td>au</td></tr><tr><td>e<sub>1</sub></td><td>0.0000<sup>+0.0020</sup><sub>-0.0000</sub></td><td>e<sub>2</sub></td><td>0.1974<sup>+0.0009</sup><sub>-0.0009</sub></td><td>e<sub>3</sub></td><td>0.569<sup>+0.022</sup><sub>-0.023</sub></td><td></td></tr><tr><td>i<sub>1</sub></td><td>86.5<sup>+3.8</sup><sub>-1.5</sub></td><td>i<sub>2</sub></td><td>86.7<sup>+2.1</sup><sub>-1.9</sub></td><td>i<sub>3</sub></td><td>-26.3<sup>+11.6</sup><sub>-8.0</sub></td><td>deg</td></tr><tr><td>Ω<sub>1</sub></td><td>331.4<sup>+1.4</sup><sub>-2.0</sub></td><td>Ω<sub>2</sub><sup>†</sup></td><td>328.9<sup>+1.4</sup><sub>-1.2</sub></td><td>Ω<sub>3</sub></td><td>108.3<sup>+3.3</sup><sub>-3.3</sub></td><td>deg</td></tr><tr><td>ω<sub>1</sub></td><td>274.11<sup>+0.15</sup><sub>-0.15</sub></td><td>ω<sub>2</sub></td><td>9.62<sup>+0.14</sup><sub>-0.14</sub></td><td>ω<sub>3</sub></td><td>9.0<sup>+3.2</sup><sub>-3.2</sub></td><td>deg</td></tr><tr><td>M<sub>1</sub></td><td>176.02<sup>+0.15</sup><sub>-0.15</sub></td><td>M<sub>2</sub></td><td>85.68<sup>+0.13</sup><sub>-0.12</sub></td><td>M<sub>3</sub></td><td>31.3<sup>+1.4</sup><sub>-1.4</sub></td><td>deg</td></tr><tr><td>γ</td><td>8.5<sup>+1.6</sup><sub>-1.6</sub></td><td></td><td></td><td></td><td></td><td>km s<sup>-1</sup></td></tr><tr><td>i'<sub>1</sub></td><td>93.5<sup>+1.5</sup><sub>-2.8</sub></td><td></td><td></td><td></td><td></td><td>deg</td></tr><tr><td>Ω'<sub>1</sub></td><td>148.5<sup>+5.8</sup><sub>-2.9</sub></td><td></td><td></td><td></td><td></td><td>deg</td></tr><tr><td></td><td></td><td></td><td></td><td>i'<sub>3</sub></td><td>25.6<sup>+9.3</sup><sub>-15.2</sub></td><td>deg</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** The epoch is $T_0 = 2.456224.724705$. The values below correspond to that of best-fit solutions with $\chi^2 = 2578$, with individual contributions: $\chi^2_{rv} = 2237$, $\chi^2_{etv} = 151$, $\chi^2_{evd} = 3.3$ and $\chi^2_{sky} = 185$. The masses of components in $M_\odot$ units were $m^{Aa} = 2.232911^{+0.000091}_{-0.000093}$, $m^{Ab} = 2.009948^{+0.000092}_{-0.000093}$, $m^B = 3.7343^{+0.0070}_{-0.0073}$ and $m^C = 0.90^{+0.73}_{-0.04}$. Component B is on its lower limit $m_{\min}^B$. The mass of component C is very poorly constrained, it is more of a distant test mass compared to the others. The $3\sigma$ uncertainties of the elements were determined by a simplified one-dimensional $\chi^2$ mapping, assuming a relative increase of $\chi^2$ by a factor of 1.13, i.e. suitable for the number of degrees of freedom we have ($\nu = 908$). The uncertainty values were verified using the bootstrap method with 100 random selections of datasets and corresponding simplex optimisation, but realistic uncertainties are likely to be larger than that because there are a number of local minima with statistically equivalent $\chi^2$ values. We do not report a full correlation matrix of our solution here. Its non-diagonal terms indicate higher values of uncertainties for those elements that are strongly correlated or anti-correlated with others (e.g. $r_{m^{Aa},\omega_2} = 0.74$, $r_{a_1,\Omega_2} = -0.77$, $r_{a_1,\omega_2} = -0.80$, $r_{a_3,i_3} = -0.79}$). (f) The value is expressed in hierarchical Jacobian elements, with respect to Aa+Ab barycentre because this pair is the most compact and massive. If the reference point were the photocentre of the brightest component B instead, then the longitude of the ascending node would be shifted by $-180^\circ$.

We clearly see that the value of $\chi^2 = 2578$ is still about three times higher than the number of degrees of freedom, $\nu = N_{data} - M_{free} = 931 - 23 = 908$, and formally speaking, we should be ready to admit that our model is plainly wrong. Nevertheless, the residua seem to be distributed normally, and realistic uncertainties (including some systematics) may be larger than expected. To obtain $\chi^2 \approx \nu$ we would need measurement uncertainties as large as $\sigma_{rv} \approx 3.5$ km s$^{-1}$, $\sigma_{etv} \approx 10$ min, $\sigma_{sky} \approx 1$ mas (for component B) or 10 mas (for component C). We consider these numbers to be quite realistic given the heterogeneous data set we have. Additional problems may contribute to the error budget, such as nightly and night-to-night variations of dispersion relations, unaccounted blending of spectral lines, systematics due to the normalisation procedure, or photocentre motions of the inner binary affecting astrometric positions.

## 8.3. Differences between traditional and N-body models

Most importantly, orbital elements do change in the course of time; especially $i_1, \Omega_1, \omega_1, \Omega_2, \omega_2$ seem to be critical in the case of $\xi$ Tau (see Fig. 11). While the precession of $\omega_2$ was accounted for, the remaining terms were not. The precession of nodes $\Omega_1, \Omega_2$ about the total angular momentum axis occurs with a $\approx 19$ yr period. In the Laplace plane, which is perpendicular to the total angular momentum, this would cause a circulation of $\Omega$'s from $0^\circ$ to $360^\circ$, but we can only see an oscillation of at most $3.5^\circ$ that is due to the purely geometrical projection to the plane of the sky. There are also inevitable coupled oscillations of inclinations, with $i_1$ ranging from $84.5^\circ$ to $88.2^\circ$. All these rather expected secular effects are discussed in much more detail in Sect. 9.1.
<figure><img src="image_12.png" /></figure>

**Fig. 11.** Time evolution of the osculating orbital elements over a time span −11 000 to +1000 days from the epoch T₀ = 2 456 224.724705, covered by observations of ξ Tau. **Left:** the semi-major axis *a*₁, eccentricity *e*₁, inclination *i*₁, longitude of ascending node Ω₁, and the argument of pericentre ω₁ (poorly defined because *e*₁ → 0) of the inner, eclipsing binary orbit (components Aa and Ab). **Right:** the same parameters *a*₂, *e*₂, *i*₂, Ω₂, ω₂ for orbit 2 (i.e. components (Aa+Ab) and B). All these plots correspond to the simulation with χ² = 2578, presented in Fig. 9. Variations in the inclination *i*₁ and argument of pericentre ω₂ are of major interest, since they result in observable effects. On the other hand, the distant orbit 3 (not shown here) exhibits only minor variations of its elements. The bump in the osculation elements of orbit 2 at JD ≈ 2 455 500 is related to the passage of component C through its pericentre.

Additionally, there are short-period oscillations not de-
scribed by the secular theory. While *a*₁ and *a*₂ only oscillate
about constant mean values, there seems to be a mid-term evolu-
tion of both *e*₁ and *e*₂, with amplitudes reaching 0.008, which
is larger than the uncertainty of their initial values, that is,
*e*₂ = 0.1974<sup>+0.0009</sup><sub>-0.0010</sub>. In this particular case, this is related to the
periastron passage of component C.

We emphasise that it is absolutely necessary to use an *N*-body model (like ours), otherwise traditional methods assuming constant orbital elements (or precessing *ω*s only) may result in systematic discrepancies or artefacts. When the parameters reported in Table 15 are compared to those derived by classical models (Table 12), the general agreement between the elements is evident, but their uncertainty intervals do not always overlap. This is probably to be expected because we compare osculating (apples) and fixed orbital elements (oranges).

An outstanding example of how classical methods may fail
is a detailed analysis of MOST light curves and the correspond-
ing minima timings from 2012. At first, we thought that the un-
even spacing of minima indicates a non-zero eccentricity of the
inner orbit, *e*₁ ≈ 0.002. However, this is in stark contrast with
past RV measurements, which constrain forcing of *e*₁(*t*) due to

perturbations by component B and require *e*₁(*t* = *T*₀) → 0.
Figure 12 shows upon close scrutiny that the oscillation of the
semi-major axis *a*₁ has a period 3.76 days, which is half of the
synodic period *P*<sub>syn1</sub> of orbit 1, in a system that corotates with
orbit 2. Moreover, its amplitude slightly decreases as compo-
nent B moves farther away. These tiny perturbations are the real
cause of the observed eclipse timing variations. They also allow
us to discard mirror models with *Ω*'<sub>1</sub> ≠ *Ω*₂ and prefer those with
*Ω*₁ ≈ *Ω*₂ because the resulting χ²<sub>ev</sub> = 390 vs. 150 is significantly
different. Again, the eclipse variations are explained in more de-
tail in Sect. 9.2.

8.4. *Model with closure phases to resolve mirror solutions*

The admissible solutions presented in Table 15 are degenerate
in the sense that we cannot distinguish among several mirror
models (in particular *i*'<sub>1</sub>, *i*'<sub>3</sub>). To resolve this degeneracy, we con-
structed an *N*-body model that accounts for interferometric vis-
ibilities and closure phases. The latter are especially suitable to
detect any asymmetries, while the former are necessary to cor-
rectly obtain (symmetric) angular positions and separations.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparison of the osculating semi-major axis $a_1$ (bottom) and eccentricity $e_1$ evolution (middle) as computed by our N-body model for two mirror solutions with $\Omega_1 \approx 331^\circ$ (bold solid) and $\Omega'_1 \approx 151^\circ$ (red dashed). Only a short time span of 12 days is shown, close to the epoch $T_0$. The corresponding ETVs of minima observed by MOST are also shown at the top. The former solution $\Omega_1 \approx 331^\circ$ has the corresponding $\chi^2_{\text{etv}}$ (for all $N_{\text{etv}} = 35$ measurements) significantly lower than the latter, 150 vs. 390, so that we consider it as the preferred value.</figcaption></figure>

In addition to Eqs. (14) to (18), we have a few more relations:

$$V'(u, v) = \frac{1}{L_{\text{tot}}} \sum_{j=1}^{4} L_{ij} 2 \frac{J_1(\pi\theta_j\sqrt{u^2+v^2})}{\pi\theta_j\sqrt{u^2+v^2}} e^{-2\pi i(ux'_{aj}+vy'_{aj})}, \quad (21)$$

$$L_{ij}(T_{\text{eff}j}, R_j) \approx \int_{\lambda_i - \Delta\lambda_i/2}^{\lambda_i + \Delta\lambda_i/2} 4\pi R_j^2 \pi B_{\lambda}(T_{\text{eff}j}) \, d\lambda, \quad (22)$$

$$\chi^2_{\text{vis}} = \sum_{i=1}^{N_{\text{vis}}} \frac{(|V'(u_i, v_i)|^2 - |V_i|^2)^2}{\sigma^2_{\text{vis}i}}, \quad (23)$$

$$T'_3 = V'(u_1, v_1)V'(u_2, v_2)V'(-u_1 + u_2), -v_1 + v_2), \quad (24)$$

$$\chi^2_{\text{clo}} = \sum_{i=1}^{N_{\text{clo}}} \frac{(\arg T'_3 - \arg T_{3i})^2}{\sigma^2_{\text{clo}i}}, \quad (25)$$

with the notation described in Table 16. The complex visibilities $V'$ and their triple products $T'_3$ were computed assuming uniform disks for individual components. Relative luminosities $L_{ij}$ at a given effective wavelength $\lambda$ were computed by a black-body approximation.

This extended model minimises $\chi^2 = \chi^2_{\text{rv}} + \chi^2_{\text{etv}} + \chi^2_{\text{edv}} + \chi^2_{\text{sky}} + \chi^2_{\text{vis}} + \chi^2_{\text{clo}}$ and has nine additional free parameters: distance $d$ to $\xi$ Tau, uniform-disk radii $R_j$, and effective temperatures $T_{\text{eff}j}$ of all the components, even though the contribution of component C is only minor (clearly lower than 10% at the longest wavelength, $\lambda = 2.6 \, \mu\text{m}$).

We used all observational data from the MARKIII, NPOI, CHARA/VEGA, and VLTI/AMBER spectro-interferometers,

<table><caption>**Table 16.** Notation used for additional coordinates and quantities needed in our extended N-body model.</caption><tr><td>x<sub>a</sub>, y<sub>a</sub></td><td>Aa-centric angular coordinates</td></tr><tr><td>V</td><td>complex visibility; squared visibility is |V|^2</td></tr><tr><td>T<sub>3</sub></td><td>complex triple product; closure phase is arg T<sub>3</sub></td></tr><tr><td>u, v</td><td>projected baselines (expressed in cycles, B/λ)</td></tr><tr><td>θ = 2R/d</td><td>angular diameter</td></tr><tr><td>d</td><td>distance to the system</td></tr><tr><td>L, L<sub>tot</sub></td><td>component luminosity and the total luminosity</td></tr><tr><td>T<sub>eff</sub></td><td>effective temperature</td></tr><tr><td>R</td><td>stellar radius (uniform disk)</td></tr><tr><td>λ, Δλ</td><td>effective wavelength and bandwidth</td></tr><tr><td>B<sub>λ</sub>(T)</td><td>the Planck function</td></tr><tr><td>σ<sub>vis</sub></td><td>uncertainty of the squared visibility</td></tr><tr><td>σ<sub>clo</sub></td><td>uncertainty of the closure phase</td></tr></table>

with $N_{\text{vis}} = 17391$ measurements of the squared visibility $|V|^2$ and $N_{\text{clo}} = 4856$ measurements of the closure phase $\arg T_3$ (from NPOI and VLTI/AMBER). The total number of degrees of freedom is thus $\nu = N_{\text{data}} - M_{\text{free}} = 28019 - 32 = 27987$. At the same time, we did not use astrometric positions ($\chi^2_{\text{sky}}$) of component B because they are not independent; all the information should be contained in $|V|^2$ and $\arg T_3$ measurements.

Initially, we used nominal uncertainties and weights $w_{\text{vis}} = 1$, $w_{\text{clo}} = 1$, but the resulting $\chi^2_{\text{vis}} + \chi^2_{\text{clo}}$ value was too high ($\approx 10^5$), even for our best-fit models (cf. Fig. 13). The most likely reason is that we did not account properly for all calibration uncertainties. To resolve this problem, an internal re-calibration would be necessary. A possible explanation for the too high $\chi^2$ has been given in Sect. 6. For example, CHARA/VEGA interferometry from Sept. 29 2012 exhibits unrealistically quick changes of $|V|^2$ at an almost constant baseline $B/\lambda \approx 1.3$ to $1.4 \times 10^8$ cycles (see Fig. C.8). In our case, we decreased the weight $w_{\text{vis}} = 0.1$ to avoid it dominating other $\chi^2$ contributions (e.g. eclipse timing variations).

We focused on a limited set of seven mirror models, always with one or two modified orbital elements (see Table 17). For each of them, we performed one simplex run, verified by simulated annealing with the initial temperature 100 000 kelvin, schedule $T^{i+1} = 0.99T^i$ and 100 iterations at given $T^i$, so that other free parameters were able to adapt themselves to a new situation, and we computed $\chi^2$s that are reported in the same table. If the final value remains relatively high, it means the model is not compatible with the respective interferometric data.

Clearly, we are sufficiently sensitive to resolve $\Omega_2$ and $i_2$, that is, the longitude of the ascending node and the inclination of component B (see Fig. 14), but not directly to resolve $\Omega_1$, $i_1$, or $i_3$ elements. Consequently, we can discard $\Omega'_2$, $i'_2$ and prefer $\Omega_2 \approx 331^\circ$, $i_2 \approx 86^\circ$ solution on the basis of the closure phase measurements alone.

Moreover, because our N-body model is constantly constrained by RV, ETV, ETD, and astrometric data, which prevent a convergence to unrealistic values of all the parameters, we can spot (in Table 17) that the squared visibility measurements are not compatible with $\Omega'_1$ and $i'_1$, therefore they were discarded as well and the $\Omega_1 \approx 329^\circ$, $i_1 \approx 86^\circ$ solution was preferred.

Finally, as demonstrated in Sect. 8.3, the N-body dynamics and ETV measurements allow us to safely discard any $\Omega_1 \neq \Omega_2$, therefore we clearly prefer $\Omega_1 = 329^\circ$. The only remaining ambiguity is thus the inclination $i_3$ vs. $i'_3$. To conclude this section, a combination of approximately orthogonal measurements (RV, ETV, ETD, $|V|^2$, $\arg T_3$) leads to interesting and solid results, which is as expected.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Distributions of normalised residuals ($|V'|^2 - |V|^2/\sigma_{vis,i}$) of the squared visibility for our best-fit model with $\chi^2_{vis} = 133492$, while the total number of measurements is $N_{vis} = 17391$. Four separate datasets are shown, corresponding to the MARKIII, NPOI, CHARA/VEGA, and VLTI/AMBER interferometers. The distributions are not perfectly symmetric about zero and for CHARA/VEGA data are significantly wider, probably due to unaccounted calibration uncertainties.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Distributions of normalised residuals ($\arg T'_3 - \arg T_{3i})/\sigma_{clo,i}$) of the closure phase measured with the NPOI instrument for two best-fit models with different values of the longitude of the ascending node $\Omega_2 = 329^\circ$ and $\Omega'_2 = 149^\circ$. Both distributions seem symmetric about the origin, indicating there are no serious systematics in $\arg T_3$ measurements. However, the former distribution is substantially narrower than the latter, so that the mirror solution $\Omega'_2$ can be discarded.</figcaption></figure>

<table><caption>**Table 17.** Summary of $\chi^2_{vis}$ and $\chi^2_{clo}$ values for squared visibility $|V|^2$ and closure phase $\arg T_3$ measurements.</caption><thead><tr><th>Orbital elements</th><th>χ<sup>2</sup><sub>vis</sub></th><th>χ<sup>2</sup><sub>clo</sub></th><th>Note</th></tr></thead><tbody><tr><td>nominal</td><td>133492</td><td>23151</td><td></td></tr><tr><td>Ω'<sub>1</sub> = Ω<sub>1</sub> + 180° ≈ 151°</td><td><b>162 632</b></td><td>23053</td><td>×</td></tr><tr><td>Ω'<sub>2</sub> = Ω<sub>2</sub> + 180° ≈ 149°</td><td>355456</td><td><b>105 975</b></td><td>×</td></tr><tr><td>(Ω'<sub>1</sub>, Ω'<sub>2</sub>)</td><td>322079</td><td><b>100 480</b></td><td>×</td></tr><tr><td>i'<sub>1</sub> = 180° - i<sub>1</sub> ≈ 94°</td><td><b>149 901</b></td><td>24683</td><td>×</td></tr><tr><td>i'<sub>2</sub> = 180° - i<sub>2</sub> ≈ 94°</td><td>734267</td><td><b>69 102</b></td><td>×</td></tr><tr><td>i'<sub>3</sub> = -i<sub>3</sub> ≈ 24°</td><td>138316</td><td>23393</td><td></td></tr><tr><td>(i'<sub>1</sub>, i'<sub>2</sub>)</td><td>755013</td><td><b>69 463</b></td><td>×</td></tr></tbody></table>

**Notes.** Only a limited set of mirror models is shown with respect to the nominal model ($\Omega_1 \approx 331^\circ$, $\Omega_2 \approx 329^\circ$, $\Omega_3 \approx 110^\circ$, $i_1 \approx 86^\circ$, $i_2 \approx 86^\circ$, $i_3 \approx -24^\circ$). The closure phase measurements allow us to discard four of them, namely those with $\Omega'_2$ and $i'_2$, because the $3\sigma$ level corresponds to a relative increase by 1.051, i.e. $\chi^2_{clo} \approx 24331$. Moreover, the $|V|^2$ measurements do not favour $\Omega'_1$ and $i'_1$ ($3\sigma$ is at 1.028, $\chi^2_{vis} \approx 137229$). The symbol $\times$ in the last column indicates we discard this possibility.

We also comment on the fact that even this type of model may be insufficient. Other physical effects exist that we did not account for, such as tidal interactions of non-spherical stars, spin-orbital coupling, various magneto-hydrodynamic phenomena, or pulsations of (all) components. Their importance for the dynamics of ξ Tau is yet to be assessed.

## 9. Dynamical evolution of the Aa+Ab+B subsystem

The osculating orbital elements shown in Fig. 11 exhibit many variations over different timescales, from the short period of the inner eclipsing binary, to the intermediate period of the orbital motion of component B with respect to the eclipsing binary, up to long periods of tens to hundreds of years. Are we able to understand some of these terms, including their amplitude, and determine parametric dependencies on stellar masses and periods of orbits 1 and 2? To do so, we need to turn to perturbation theory. In this section we neglect dynamical effects of the distant component C and focus on the triple subsystem Aa+Ab+B.

The hierarchy of the ξ Tau system implies a preferential choice of Jacobi coordinates to describe its dynamics, in which

on the one hand, *r* is the relative position of Ab with respect to Aa, and on the other, *R* is the relative position of component B with respect to the barycentre of orbit 1. The conjugate momenta involve reduced masses $m'_1 = m^{Aa}m^{Ab}/M^1$ and $m'_2 = m^{B}M_1/M_2$ of orbits 1 and 2, with $M_1 = m^{Aa} + m^{Ab}$ and $M_2 = M_1 + m^{B}$. To zero-order approximation, both systems evolve on Keplerian orbits, but their interaction introduces a perturbation that causes *r* and *R* to follow trajectories described by numerical integrations in Sect. 8 that are more complicated than Keplerian orbits. The elliptical approximation may be only applicable to a certain interval of time. The latter becomes short especially for compact systems. ξ Tau is a good representative of this class.

In the world of perturbation theory, both orbits 1 and 2 are represented by a set of osculating orbital elements that evolve in time as a result of their mutual interaction. From a plethora of perturbations described in this way, we recall two results relevant for the observed features of the ξ Tau system. We first describe the secular effects, whose duration is conveniently short for this compact system to be detected, and then some of the long- and short-period eclipse time variations in the eclipsing binary.

### 9.1. Secular effects

We define Delaunay momenta $L_1 = m'_1 \sqrt{GM_1a_1} = m'_1 n_1 a_1^2$ and $L_2 = m'_2 \sqrt{GM_2a_2} = m'_2 n_2 a_2^2$ of orbits 1 and 2 (e.g. Harrington 1968, 1969; Soderhjelm 1975; Breiter & Vokrouhlický 2015). Here $n_1$ and $n_2$ are the mean motion values of the orbits 1 and 2, both related to the semi-major axes $a_1$ and $a_2$ through the third Kepler law: $n_1^2 a_1^3 = GM_1$ and $n_2^2 a_2^3 = GM_2$ ($G$ is the gravitational constant). In a secular approximation, when the orbital longitude for both orbits 1 and 2 is removed from the interaction (e.g. Harrington 1969; Breiter & Vokrouhlický 2015), the semi-major axes $a_1$ and $a_2$ are constant.

The dynamics of the Aa+Ab+B system may in principle be studied in an arbitrary reference frame. However, its description becomes very simple in a preferred frame that is often called Laplacian. The z-axis of this frame is aligned with the total orbital angular momentum of the system. To distinguish osculating orbital elements in the observer-oriented frame we used above, we denote the elements in the Laplacian frame with a tilde. For instance, the orbital inclinations for orbits 1 and 2 are denoted $\tilde{i}_1$ and $\tilde{i}_2$, and the corresponding longitudes of nodes $\tilde{\Omega}_1$ and $\tilde{\Omega}_2$.
The secular evolution of the triple system is particularly simple when three conditions are met: the eccentricity $e_1$ of the inner orbit is negligible, the mutual angle $\tilde{J} = \tilde{\tau}_1 + \tilde{\tau}_2$, of orbital planes 1 and 2 is small, and the system is wide enough, such that on the timescale of interest only the quadrupole interaction of the inner and outer orbits is relevant. The mutual angle $\tilde{J}$ can be determined by the orbital elements in the observer reference frame using

$$ \cos \tilde{J} = \cos i_2 \cos i_1 + \sin i_2 \sin i_1 \cos (\Omega_1 - \Omega_2). \quad (26) $$

These conditions fortunately currently apply to the $\xi$ Tau system. We also note for the third condition that the octupole interaction is very small because of nearly equal masses in orbit 1, i.e., $m^{\text{Aa}} \approx m^{\text{Ab}}$. The next secular contribution would arise from the non-linear quadrupole effect (e.g. Breiter & Vokrouhlický 2015), which is small on a timescale of some decades. Then, $e_1 = 0$ is a stable solution, and $e_2$ and $\tilde{J}$ are constant in time. When the orbital elements are referred to the invariable plane that is normal to the total angular momentum, the orbital inclinations $\tilde{\tau}_1$ and $\tilde{\tau}_2$ of orbits 1 and 2 are constant as well, and both orbital planes uniformly precess in the inertial space about the total angular momentum direction. Their nodes $\tilde{\Omega}_1 = \tilde{\Omega}$ and $\tilde{\Omega}_2 = \tilde{\Omega} + \pi$ linearly advance with a rate (e.g. Soderhjelm 1975; Breiter & Vokrouhlický 2015)

$$ \frac{\dot{\tilde{\Omega}}}{n_2} \approx \frac{3}{4\eta_2^3} \frac{m^B}{M_2} \frac{n_2}{n_1} \cos \tilde{J} \sqrt{1 + \gamma^2 + 2\gamma \cos \tilde{J}}, \quad (27) $$

where $\gamma = L_1/(L_2\eta_2)$ is the ratio of the angular momenta of the two orbits, and $\eta_2 = \sqrt{1 - e_2^2}$. In triple systems the outer orbit typically has a dominant share of the total angular momentum of the system, thus $\gamma < 1$. For $\xi$ Tau we have approximately $\gamma \approx 0.132$. Unless precisely coplanar, the main effect of the orbital-plane precession is in periodic changes of inclinations $i_1$ and $i_2$ in the observer system. These variations directly affect magnitude depths of the eclipses, or might eventually cause the system to become non-eclipsing for a certain period of time.

In addition to the steady precession of the orbital planes, the second secular effect in the given setup consists of precession of the pericentre of the outer orbit. Denoting its longitude $\tilde{\omega}_2$, we have

$$ \frac{\dot{\tilde{\omega}}_2}{n_2} \approx \frac{3}{8\eta_2^3} \frac{m^B}{M_2} \frac{n_2}{n_1} \gamma \left( 3 \cos^2 \tilde{J} - 1 - \frac{\gamma \sin \tilde{J} \sin 2\tilde{J}}{1 + \gamma \cos \tilde{J} + \sqrt{1 + \gamma^2 + 2\gamma \cos \tilde{J}}} \right). \quad (28) $$

Comparing Eqs. (27) and (28), we note that the pericentre precession frequency of the outer orbit is slower by a factor $\approx \gamma$ than the nodal frequency (assuming $\tilde{J}$ sufficiently small). Thus nodes and inclinations vary on shorter timescales than the argument of pericentre of orbit 2.

## 9.2. Long- and short-period eclipse variations

The mutual interaction of the orbits also results in a palette of periodic perturbations. So far, the long-period effects, namely those with a period $P_2$ of orbit 2, have been extensively studied (e.g. Soderhjelm 1975, 1982; Borkovits et al. 2003, 2011, 2015). We here focus on the ETVs, that is, on advances and delays $\delta t_{LP}$ in epochs of eclipse of orbit 1 that are due to the variations in

its mean motion $n_1$ caused by component B. Assuming for simplicity coplanar orbits $\tilde{J} = 0$ deg, we obtain (e.g. Soderhjelm 1975; Borkovits et al. 2011, 2015; Rappaport et al. 2013, which also contain terms proportional to $\propto \tilde{J}^2$)

$$ n_1 \delta t_{LP} \approx \frac{m^B}{M_2} \frac{n_2}{n_1} W(e_2, \ell_2), \quad (29) $$

with

$$ W(e_2, \ell_2) = f_2 - \ell_2 + e_2 \sin f_2, \quad (30) $$

where $f_2$ and $\ell_2$ are the true and mean anomalies of orbit 2. For a low eccentricity $e_2$ we have $W \approx 3e_2 \sin \ell_2$. Obviously, the principal component of ETV in Eq. (29) becomes zero for a circular orbit 2 because it is related to variations of $n_1$ triggered by variations in the distance $R$ to component B.

In the course of this work, we noted that dominant short-period effects may also be of interest (those with the period of the inner orbit 1), provided high-quality eclipse data are collected. Using methods of first-order perturbation theory, we found that the leading short-period term reads

$$ n_1 \delta t_{SP} \approx \frac{21}{8} \frac{m^B}{M_2} \left( \frac{n_2}{n_1} \right)^2 \left( \frac{a_2}{R} \right)^3 \sin 2(\lambda_1 - \tilde{F}_2), \quad (31) $$

where $R$ is the distance of component B to the barycentre of the inner binary system, and $\tilde{F}_2 = \tilde{\omega}_2 + f_2$ is its true orbital longitude. The term has a period equal to half the synodic period of the Aa+Ab system in a reference frame corotating with the motion of component B.

This effect is not primarily dependent of the eccentricity $e_2$ because it is triggered by variations in the mutual positions of components Aa and Ab with respect to component B. Its magnitude is smaller by a factor 0.4 at periastron and by 0.1 at apoastron of orbit 2. Nevertheless, the effect is not entirely negligible, and we found that it contributes to the observed eclipse shift in the MOST data (see Fig. 12).

## 9.3. Comparison of the secular theory with the results of the analytic and numerical models

Here we compare the apsidal motion detected with both analytic and numerical methods and additional secular and periodical variations of orbital elements predicted by the numerical model presented in Sect. 8

- The apsidal motion of orbit 2: first, we use results of the analytic theory above. Using nominal orbital parameters from Table 12, we obtain $\tilde{J} = 0.19 \pm 1.89$ deg, and consequently $\dot{\omega}_2 = 2.185 \pm 0.045$ deg yr$^{-1}$. We note that $\dot{\omega}_2$ may be directly obtained from Eq. (28) because the nodal longitude $\tilde{\Omega}_2$ in the observer frame oscillates without any secular drift. This is about a third lower than the value detected with the analytic RV curve model (see Table 5), but in excellent agreement with the N-body model, whose prediction is $\dot{\omega}_2 = 2.11$ deg yr$^{-1}$, and with fit of the interferometric observations (see Table 11).

- The nodal motion of orbits 1 and 2: inserting nominal parameters from Table 12 provides the mean nodal drift $\tilde{\Omega} = 18.98 \pm 0.53$ deg yr$^{-1}$ (Eq. (27)), which is again in excellent agreement with results of the N-body model; we note that the periods of the nodal oscillations are effectively $\approx 19.43$ deg yr$^{-1}$ for orbit of component A ($\Omega_1$) and $\approx 19.81$ deg yr$^{-1}$ for orbit B ($\Omega_2$). Values are not exactly the
same, probably because of interaction with component C, which was not included in the perturbation theory. There is a hint of a shallower depth of the Hvar photometric observations from early 2007 when our model predicts a higher value of the inclination $i_1$. However, to determine the inclination variations, we need more accurate photometric observations in the future.

- *Eclipse-timing variations – orbit 1:* Eqs. (29)–(31) provide amplitudes of the ETVs (assuming that component B is at periastron) of orbit 1 $\delta t_{ETV, long} = 0.0162 \pm 0.0007$ d, $\delta t_{ETV, short} = 0.0068 \pm 0.0003$ d. Their sum agrees with the detected amplitude of ETVs ($\delta t_{ETV, OBS} = 0.025 \pm 0.010$ d). We also note that the two primary eclipse minima in the MOST data were found to be phase-shifted by ≈0.0003 in Sect. 4.3. This is about 0.1° in orbital longitude of inner orbit 1. Combining results in Eqs. (29)–(31) and taking into account $\ell_2 \approx 86^\circ$ and $\lambda_1 \approx F_2$ from Table 15, we obtain very good agreement with the observed shift.

# 10. Motivations for future observations of ξ Tauri

First, it seems desirable to continue the observations of the times of minima and, more importantly, eclipse duration and depth. At an epoch after approximately RJD 59 405.0, that is, in the second half of 2021, we would expect either persisting or disappearing eclipses of the inner pair Aa+Ab for different mirror solutions. Consequently, this is a direct and independent test of our analysis of closure phase measurements in Sect. 8.4. We note that the nominal solution shown in Fig. 11 exhibits too small variations of $i_1$, such that the eclipsing binary would be eclipsing constantly.

Nevertheless, even the nominal solution predicts nearly 4° full amplitude of variation in $i_1$ and we expect fairly well observable effects. We suggest, for instance, a space-born observation of a similar quality to that of MOST, obtained at the turn of 2016 and 2017, when the predicted $i_1$ value would be highest (about 88.2°). The change in eclipse depth, as compared to the MOST data, should be about 0.05 mag, which is very easily detectable. Such a single observation would further constrain parameters of ξ Tau with an exceptional accuracy.

It would be of great help if the line spectra of the faint component C, separated by 200 to 600 mas from the triple Aa+Ab+B, were obtained and the corresponding RV measured. This would also allow us to distinguish between the remaining two mirror orbital solutions for the motion of this component.

Precise and uninterrupted space-based photometry on a longer time-span would be useful to unambiguously resolve oscillation modes and splittings. Given the high rotation frequency $f_{rot} = 2.38$ d$^{-1} \approx 27.5$ µHz, it should not be that difficult (the minimum time-span $\Delta t \approx 1/f_{rot}$), but currently aliases with instrumental frequencies seem to limit the S/N in the Fourier spectrum.

As an alternative, series of high-resolution high S/N spectra would be needed to detect the oscillation modes independently, as the travelling sub-features in the line profiles of component B are broadened by a relatively high rotation. Precise RV measurements of components Aa and Ab may also reveal the Rossiter-McLaughlin effect, which gives the rotational sense of the two components.

A new series of long-baseline optical spectro-interferometric observations including measurements of closure phase are highly desirable, because they would provide a fully independent estimate of the orbital elements of orbit 1, would independently determine the sense of revolution of the components of orbit 1

with respect to orbit 2, and would finally provide an independent estimate of the radii of components Aa and Ab.

# 11. Conclusion

We have conducted an in-depth study of the quadruple stellar system ξ Tau, starting from simple analytic models for different types of observations (see Sects. 3–7), and concluding with a complex *N*-body model that combines astrometric, photometric, spectroscopic, and spectro-interferometric observations to a certain degree (see Sect. 8). We were able to set tight constraints on three components of ξ Tau, and they will provide an excellent test case for models of stellar evolution, while the full description of the geometry of the hierarchy will provide a test of the binary formation.

The analytic models allowed us to estimate properties of components Aa, Ab, and B that are highly consistent (see the critical summary of the analytic models in Sect. 7) and mean orbital elements of orbits 1, 2, and 3 using different methods that are again consistent with each other, but provided limited-to-no insight into the dynamic evolution of orbits of the ξ Tau.

This discrepancy was fixed with the *N*-body model, which properly accounts for the dynamic interaction within the system and is able to fit RVs, ETVs, and astrometric positions simultaneously. It provided a set of osculating elements and component masses whose evolution fits the observables (see Table 13). It also provided insight into the long- and short-term evolution of the osculating elements (see Fig. 11) and also resolved thegrade and retrograde solution (between orbits 1 and 2) solely from ETVs. The result also supports previous analyses because it did not vary much from their outcome.

Perturbation theory shows that the most pronounced secular evolution of elements, that is, the advance of the apsidal line of orbit 2, the harmonic variation of the inclination $i_{1,2}$, and the longitude of the ascending node $\Omega_{1,2}$, are explained by a quadrupole interaction between orbits 1 and 2. The same applies to the predicted size of ETVs, which agree well with observations.

*Acknowledgements.* The research of J.N., P.H., M.W., and P.Z. was supported by grants P209/10/0715 and GA15-02112S of the Czech Science Foundation. The research of J.N. and P.H. was also supported by grant No. 678212 of the Grant Agency of the Charles University. The Navy Prototype Optical Interferometer is a joint project of the Naval Research Laboratory and the US Naval Observatory, in cooperation with Lowell Observatory and is funded by the Office of Naval Research and the Oceanographer of the Navy. The authors thank Jim Benson and the NPOI observational support staff, whose efforts made this project possible. This research has made use of the SIMBAD astronomical literature database, operated at CDS, Strasbourg, France. The CHARA Array is operated with support from the National Science Foundation through grant AST-0908253, the W. M. Keck Foundation, the NASA Exoplanet Science Institute, and from Georgia State University. This publication is supported as a project of the Nordrhein-Westfälische Akademie der Wissenschaften und der Künste in the framework of the academy programme by the Federal Republic of Germany and the state Nordrhein-Westfalen. H.B. acknowledges financial support by Croatian Science Foundation under the project 6212 "Solar and Stellar Variability". The project is based on data obtained from the ESO Science Archive Facility under request number jnemravova217453, on spectral data retrieved from the ELODIE archive at Observatoire de Haute-Provence (OHP), and on observations made at the South African Astronomical Observatory (SAAO). P.Z. wish to thank the staff at SAAO for their warm hospitality and help with the equipment. A.F.J.M. is grateful for financial assistance to NSERC (Canada) and FQRNT (Quebec). The observations with the MPG 2.2 m telescope were supported by the Czech Ministry of Education, Youth and Sports project LG14013 ("Tycho Brahe: Supporting Ground-based Astronomical Observations") during run P2 in May 2015. We acknowledge the use of the electronic database from the CDS, Strasbourg and electronic bibliography maintained by the NASA/ADS system. We acknowledge the constructive criticism by the referee Peter P. Eggleton, which helped us to improve the paper.
## References

Aller, L. H., Appenzeller, I., Baschek, B., et al. 1982, Landolt-Börnstein: Numerical Data and Functional Relationships in Science and Technology – New Series, Group 6 Astronomy and Astrophysics, Vol. 2

Armstrong, J. T., Mozurkewich, D., Rickard, L. J., et al. 1998, ApJ, 496, 550

Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

Bolton, C. T., & Grunhut, J. H. 2007, in IAU Symp. 240, eds. W. I. Hartkopf, P. Harmanec, & E. F. Guinan, 66

Bonneau, D., Clausse, J.-M., Delfosse, X., et al. 2006, A&A, 456, 789

Borkovits, T., Érdi, B., Forgács-Dajka, E., & Kovács, T. 2003, A&A, 398, 1091

Borkovits, T., Csizmadia, S., Forgács-Dajka, E., & Hegedüs, T. 2011, A&A, 528, A53

Borkovits, T., Rappaport, S., Hajdu, T., & Sztakovics, J. 2015, MNRAS, 448, 946

Breiter, S., & Vokrouhlický, D. 2015, MNRAS, 449, 1691

Brož, M., Mayer, P., Pribulla, T., et al. 2010, AJ, 139, 2258

Butterworth, S. 1930, Wireless Engineer, 7

Campbell, W. W. 1909, ApJ, 29, 224

Claret, A. 1998, A&AS, 131, 395

Claret, A. 2001, MNRAS, 327, 989

Degroote, P., Acke, B., Samadi, R., et al. 2011, A&A, 536, A82

de Laverny, P., Recio-Blanco, A., Worley, C. C., & Plez, B. 2012, A&A, 544, A126

Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A, 409, 205

Eggleton, P. P., & Kiseleva-Eggleton, L. 2001, ApJ, 562, 1012

Eggleton, P. P., & Tokovinin, A. A. 2008, MNRAS, 389, 869

ESA 1997, VizieR Online Data Catalog, I/239

Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298

Fekel, J. F. C. 1981, ApJ, 246, 879

Flower, P. J. 1996, ApJ, 469, 355

Fuhrmann, K., Chini, R., Hoffmeister, V. H., et al. 2011, MNRAS, 411, 2311

Glindemann, A., Albertsen, M., Andolfato, L., et al. 2004, in New Frontiers in Stellar Interferometry, ed. W. A. Traub, Proc. SPIE, 5491, 447

Hadrava, P. 1995, A&AS, 114, 393

Hadrava, P. 1997, A&AS, 122, 581

Hadrava, P. 2009, ArXiv e-prints [arXiv:0909.0172]

Harmanec, P. 1988, Bulletin of the Astronomical Institutes of Czechoslovakia, 39, 329

Harmanec, P. 1998, A&A, 335, 173

Harmanec, P., & Horn, J. 1998, J. Astron. Data, 4, 5

Harmanec, P., Horn, J., & Juza, K. 1994, A&AS, 104, 121

Harrington, R. S. 1968, AJ, 73, 190

Harrington, R. S. 1969, Celest. Mech., 1, 200

Hummel, C. A., Mozurkewich, D., Armstrong, J. T., et al. 1998, AJ, 116, 2536

Hummel, C. A., Benson, J. A., Hutter, D. J., et al. 2003, AJ, 125, 2630

Hummel, C. A., Zavala, R. T., & Sanborn, J. 2013, Central European Astrophysical Bulletin, 37, 127

Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, The Messenger, 95, 8

Kraft, D. 1988, A software package for sequential quadratic programming, Deutsche Forschungs- und Versuchsanstalt für Luft- und Raumfahrt Köln: Forschungsbericht (Wiss. Berichtswesen d. DFVLR)

Lafrasse, S., Mella, G., Bonneau, D., et al. 2010, VizieR Online Data Catalog, II/300

Lenz, P., & Breger, M. 2004, in The A-Star Puzzle, eds. J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss, IAU Symp., 224, 786

Levison, H. F., & Duncan, M. J. 2013, Astrophysics Source Code Library [record ascl:1303.001]

Mason, B. D., Martin, C., Hartkopf, W. I., et al. 1999, AJ, 117, 1890

Moultaka, J., Ilovaisky, S. A., Prugniel, P., & Soubiran, C. 2004, PASP, 116, 693

Mourard, D., Clausse, J. M., Marcotto, A., et al. 2009, A&A, 508, 1073

Mozurkewich, D., Armstrong, J. T., Hindsley, R. B., et al. 2003, AJ, 126, 2502

Nemravová, J. A., Harmanec, P., Bencheikh, J., et al. 2013, Central European Astrophysical Bulletin, 37, 207

Palacios, A., Gebran, M., Josselin, E., et al. 2010, A&A, 516, A13

Petrov, R. G., Malbet, F., Weigelt, G., et al. 2007, A&A, 464, 1

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1993, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd edn. (New York: Cambridge University Press)

Prša, A., & Zwitter, T. 2005, ApJ, 628, 426

Prša, A., & Zwitter, T. 2006, Ap&SS, 304, 347

Rappaport, S., Deck, K., Levine, A., et al. 2013, ApJ, 768, 33

Rica Romero, F. M. 2010, Rev. Mex. Astron. Astrofis., 46, 263

Shao, M., Colavita, M. M., Hines, B. E., et al. 1988, A&A, 193, 357

Simon, K. P., & Sturm, E. 1994, A&A, 281, 286

Soderhjelm, S. 1975, A&A, 42, 229

Soderhjelm, S. 1982, A&A, 107, 54

Steiner, I., Stahl, O., Seifert, W., Chini, R., & Quirrenbach, A. 2008, in SPIE Conf. Ser., 7014, 4

Stellingwerf, R. F. 1978, ApJ, 224, 953

Storn, R., & Price, K. 1997, Journal of Global Optimization, 11, 341

Tallon-Bosc, I., Tallon, M., Thiébaut, E., et al. 2008, in SPIE Conf. Ser., 7013, 70131J

ten Brummelaar, T. A., McAlister, H. A., Ridgway, S. T., et al. 2005, ApJ, 628, 453

Tody, D. 1986, in Instrumentation in astronomy VI, ed. D. L. Crawford, SPIE Conf. Ser., 627, 733

Tody, D. 1993, in Astronomical Data Analysis Software and Systems II, eds. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, ASP Conf. Ser., 52, 173

Tokovinin, A. A. 1997, A&AS, 124, 75

Tokovinin, A., Mason, B. D., Hartkopf, W. I., Mendez, R. A., & Horch, E. P. 2015, AJ, 150, 50

van Belle, G. T., Creech-Eakman, M. J., & Hart, A. 2009, MNRAS, 394, 1925

van Leeuwen, F. 2007, A&A, 474, 653

Walker, G., Matthews, J., Kuschnig, R., et al. 2003, PASP, 115, 1023

Zasche, P., & Wolf, M. 2007, Astron. Nachr., 328, 928

Zasche, P., Uhlář, R., Kučáková, H., Svoboda, P., & Mašek, M. 2014, IBVS, 6114, 1

<footnote n="1">¹ Astronomical Institute of the Charles University, Faculty of Mathematics and Physics, V Holešovičkách 2, 180 00 Praha 8, Troja, Czech Republic</footnote>

<footnote n="2">² Laboratoire Lagrange, OCA/UNS/CNRS UMR 7293, BP 4229, 06304 Nice Cedex, France</footnote>

<footnote n="3">³ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany</footnote>

<footnote n="4">⁴ Department of Mathematics, Physics & Geology, Cape Breton University, 1250 Grand Lake Road, Sydney, NS B1P 6L2, Canada</footnote>

<footnote n="5">⁵ Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="6">⁶ Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario, M5S 3H4, Canada</footnote>

<footnote n="7">⁷ Hvar Observatory, Faculty of Geodesy, Zagreb University, Kačićeva 26, 10000 Zagreb, Croatia</footnote>

<footnote n="8">⁸ Astronomisches Institut, Ruhr-Universität Bochum, Universitätstraße 150, 44780 Bochum, Germany</footnote>

<footnote n="9">⁹ Instituto de Astronomía, Universidad Católica del Norte, Avenida Angamos 0610, Casilla 1280 Antofagasta, Chile</footnote>

<footnote n="10">¹⁰ Department of Astronomy and Astrophysics, Villanova University, Villanova, PA 19085, USA</footnote>

<footnote n="11">¹¹ CHARA Array, Mount Wilson Observatory, Mount Wilson, CA 91023, USA</footnote>

<footnote n="12">¹² Department of Astronomy and Physics, Saint Mary's University, Halifax, N.S., B3H 3C3, Canada</footnote>

<footnote n="13">¹³ Astronomical Institute, Academy of Sciences of the Czech Republic, 251 65 Ondřejov, Czech Republic</footnote>

<footnote n="14">¹⁴ Institute of Astronomy, University Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria</footnote>

<footnote n="15">¹⁵ Department de Physique, Université de Montréal, C.P. 6128, Succursale Centre-Ville, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="16">¹⁶ Observatório do Instituto Geográfico do Exército, R. Venezuela 29 3 Esq., 1500-618, Lisbon, Portugal</footnote>

<footnote n="17">¹⁷ NASA Ames Research Center, Moffett Field, CA 94035; SETI Institute, Mountain View, CA 94043, USA</footnote>

<footnote n="18">¹⁸ Université de Lyon, Université Lyon 1, École Normale Supérieure de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR 5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="19">¹⁹ US Naval Observatory, Flagstaff Station, 10391 West Naval Observatory Road, Flagstaff, AZ 86005-8521, USA</footnote>
# Appendix A: Supplementary material to the spectroscopic observations and their analyses

Details of the reduction procedure of the spectroscopic observations used in this study along with supplementary material to its analysis are given in this section.

## A.1. Acquisition and reduction of the spectroscopic observations

The reduction procedure applied to spectra from different observatories (the labelling of observatories corresponds to that introduced in Table 2) were the following:

i) **OND**: all slit spectra were secured at the coudé focus of the 2 m reflector in Ondřejov, Czech Republic, and were recorded with the CCD detector PyLoN 2048x512BX. The bias subtraction, flat-fielding, and wavelength calibration were carried out using IRAF<fnref n="9" /> (Tody 1986, 1993). The spectra were normalised with Hermite polynomials (order $k \le 10$).

ii) **FER**: the echelle spectra were acquired with the 2.2 m ESO/MPG reflector at La Silla, Chile, and were reduced (bias subtraction, flat-fielding, and wavelength calibration) with a MIDAS pipeline developed specifically for the instrument (Kaufer et al. 1999). The studied regions of the reduced spectra were normalised with Hermite polynomials (order $k \le 10$).

iii) **BES**: the spectra were acquired with an echelle spectrograph mounted at the 1.5 m Hexapod Telescope at Cerro Amazones, Chile, which is the same as the FEROS spectrograph,

and the same MIDAS pipeline was used to carry out the reduction (bias subtraction, flat-fielding, wavelength calibration). The studied regions were normalised with Hermite polynomials (order $k \le 10$).

iv) **ELO**: the echelle spectra were obtained with the 1.93 m reflector at Observatory Haute-Provence. The initial reductions (bias subtraction, flat-fielding, wavelength calibration) were carried out with a pipeline described in Baranne et al. (1996). The studied regions were normalised with Hermite polynomials (order $k \le 10$).

v) **DDO**: the slit spectra were acquired with the 1.88 m reflector at the David Dunlap Observatory. The initial reductions (bias subtraction, flat-fielding, wavelength calibration) were carried out using IRAF. The spectra were normalised with Hermite polynomials (order $k \le 10$).

vi) **LIS**: the slit spectra were acquired with the 0.356 m reflector at the Astronomical Observatory of the Instituto Geográfico do Exército, Lisbon. The dark-frame subtraction and flat-fielding were carried out in Maxim DL<fnref n="10" />. The wavelength calibration was carried out using neon comparison spectra and telluric lines in the program Visual Spec<fnref n="11" />. The instrumental response was also removed in this program, using Castor as a reference star. The spectra were normalised with Hermite polynomials (order $k \le 10$).

## A.2. Supplementary materials to the analysis of spectroscopic observations

The spectroscopic supplementary material consists of Fig. A.1 that shows a comparison of the separated and synthetic profiles. The related analyses are described in Sects. 3.5 and 3.3.

<footnote n="9">⁹ IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.</footnote>

<footnote n="10">¹⁰ Maxim DL is a commercial software designed for astronomical imaging, http://www.cyanogen.com/maxim/main.php</footnote>

<footnote n="11">¹¹ Visual Spec is a freeware designed for the wavelength calibration and the instrumental response removal, http://www.astrosurf.com/vdesnoux/index.html</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. A.1.** Comparison of the separated and synthetic spectra. The parameters defining the synthetic spectra are listed in Table 6. In each panel we show in the top spectrum component B, in the middle spectrum component Aa, and in the bottom spectrum component Ab. The thick grey line plots spectra, the thin black line separated and re-normalised spectra, and the thin red line synthetic spectra. The residuals are computed for synthetic and re-normalised separated spectra.</figcaption></figure>
# Appendix B: Supplementary material to the photometric observations and their analysis

Details on the reduction procedure of the photometric observations used in this study along with supplementary material to their analysis are given in this section.

## B.1. Acquisition and reduction of the photometric observations

The reduction procedure applied to photometric observations from different observatories (the labelling of observatories corresponds to the labelling introduced in Table 3) were the following:

i) *HVAR*: the differential observations were obtained with the 0.65 m reflector at the Hvar Observatory, Croatia, which is equipped with a photoelectric photometer with an EMI 6256 tube. The observations were acquired relative to the comparison star 4 Tau with the check star 6 Tau observed as frequently as ξ Tau and transformed to the standard *UBV* system (*UBVR* for observations acquired after RJD = 56 000) using the non-linear transformations implemented in the reduction package HEC22¹² (see Harmanec et al. 1994; Harmanec & Horn 1998). All observations were reduced with the latest release 18, which allows for time variation in the linear extinction coefficients in the course of the observing night.

ii) *HIPP*: the all-sky observations were acquired with the 0.29 m reflector of the HIPPARCOS satellite and transformed to V magnitude using the formulae derived by Harmanec (1998).

iii) *SAAO*: the differential observations were acquired at the South African Astronomical Observatory, South Africa with 0.5 m reflector equipped with a photoelectric photometer. The observations were acquired relative to the comparison star 4 Tau and 6 Tau served as a check star and were transformed to the standard Johnson system using HEC22.

iv) *VILL*: the differential observations were acquired with the Automatic Photometric Telescope at Villanova, USA. The observations were taken relative to the comparison star 4 Tau and 6 Tau served as a check star.

v) *MOST*: the all-sky observations were obtained with the 0.15 m reflector in the MOST satellite. The initial reduction was carried out according to Walker et al. (2003) and references therein. Removal of the remaining instrumental artefacts and the stray light from Earth's atmosphere is described in Sect. 2.

## B.2. Supplementary materials to the analysis of the photometric observations

The photometric supplementary material consists of Figs. B.1 and B.2 that show the available primary and secondary light-curve minima. All minima cover a time interval no longer than 30 d. See Sect. 4 for related analyses.

<footnote n="12">¹² The whole package along with a detailed manual, auxiliary data files, and results is available at http://astro.troja.mff.cuni.cz/ftp/PHOT</footnote>
<figure><img src="image_17.png" /><figcaption>Fig. B.1. All available primary minima of orbit 1. The filters are denoted as follows: *UBV* — Johnson *UBV* filters, *MO* the MOST filter, and *A* differential measurements taken in the visible without any filter. Mean RJD is given in the bottom left corner of each panel.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. B.2.** All available secondary minima of orbit 1. The filters are denoted as follows: *UBV* – Johnson *UBV* filters, *MO* the MOST filter, A differential measurements taken in the visible without any filter. Mean RJD is given in the bottom left corner of each panel.</figcaption></figure>
# Appendix C: Supplementary material to the spectro-interferometric observations and their analyses

Details on the acquisition and reduction of the spectro-interferometric observations, along with tables and figures illustrating their analysis are presented here.

## C.1. Mark III observations

The observations were carried out using a single north-south baseline three times on January 19, October 19, and November 2, 1991. The baseline length was 32 m on the first night and 15 m on the two other nights. Visibilities were recorded in three narrow-band channels at 5000 Å, 5500 Å, and 8000 Å. μ and η Tau (limb-darkened diameters of 0.41 mas and 0.98 mas, respectively, with 10% uncertainties) served as the calibrators. The calibrated visibilities were obtained from the Mark III data archive, which was created using the reduction and calibration methods described by Mozurkewich et al. (2003).

## C.2. NPOI observations

The observations were carried out with the three-beam combiner in 1998 and 2000, and from 2003 to 2013 with the six-beam combiner. Visibilities, complex triple amplitudes, and closure phases were recorded in 16 narrow-band channels between 5500 Å and 8500 Å. The journal of the NPOI observations including the calibrator stars is given in Table C.2, and the calibrator information is given in Table C.4.

The calibrators were taken from a list of single stars maintained at NPOI with diameters estimated from V and (V-K) using the surface brightness relation by Mozurkewich et al. (2003) and van Belle et al. (2009). Values of E(B-V) were derived from comparison of the observed and theoretical colours as a function of spectral type by Schmidt-Kaler in Aller et al. (1982). Values for the extinction derived from E(B-V) were compared to estimates based on the maps by Drimmel et al. (2003), and used to correct V if they agreed within 0.5. Even though the surface brightness relation based on (V-K) colours is to first order independent of the reddening, we included this small correction. The minimum (squared) visibility amplitudes corresponding to the diameter estimates are given in Table C.4 for all NPOI observations and show that the calibrators are either unresolved or only weakly resolved.

NPOI data and their reductions were described by Hummel et al. (1998) and Hummel et al. (2003). For the first time, we used a pipeline written in GDL<fnref n="13" /> for the OYSTER<fnref n="14" /> NPOI data reduction package. The pipeline automatically edits the one-second averages produced by another pipeline directly from the raw frames, based on expected performance such as the variance of fringe tracker delay, photon count rates, and narrow-angle tracker offsets. Visibility bias corrections are derived as usual from the data recorded away from the stellar fringe packet. After averaging the data over the full length of an observation, the closure phases of the calibrators were automatically unwrapped so that their variation with time, as well as that of the visibility amplitude, could be interpolated for the observations of ξ Tau. For the calibration of the visibilities, the pipeline used all calibrator stars observed during a night to obtain smooth averages of the amplitude and phase-transfer functions using a

Gaussian kernel of 80 min in length. The residual scatter of the calibrator visibilities and phases around the average set the level of the calibration uncertainty and was added in quadrature to the intrinsic data errors. The amplitude calibration error of typically a few percent in the red channels up to 15% in the blue channels was added in quadrature to the intrinsic error of the visibilities. The phase calibration was good to about a couple of degrees.

## C.3. VEGA/CHARA observations

The observations were carried out during two runs in 2011 and in 2012. Preliminary results, based on the observations obtained during the first run, were published by Nemravová et al. (2013). The reduction procedure was the same for both runs.

Five observations were acquired in 2011. All observations were obtained in the three-telescope (3T) mode and included the CHARA baselines E1E2W2, W1W2S2, and W2E2S2, ranging from 63 m to 245 m (E1, E2, S1, S2, W1, and W2 denote the telescopes in the CHARA telescope array). Ten new observations were secured in 2012. Four of them were taken in the 3T mode and the remaining six were taken in the two-telescope (2T) mode. The 2T observations included the CHARA baselines S2S1 and E2E1, their projected lengths ranging from 34 m to 66 m. The 3T observations contained the E2E1W2 and W2W1S1 baselines, their projected lengths were from 65 m to 279 m. A detailed journal of all interferometric observations with the instrument CHARA/VEGA is in Table C.3.

The observations were obtained with two cameras centred on 5350 Å (denoted BLUE) and 7300 Å (denoted RED) at spectral resolution R ≈ 5000. Individual frames were recorded with a frequency of 100 Hz and grouped into blocks containing 2500 frames. Each block was coherently summed and each observation had from 20 to 90 blocks. Two 20 nm wide bands were chosen in the BLUE region and two 30 nm wide bands in the RED region. The four bands used are $Δλ_{IF} = \{5320-5520, 5400-5600, 7000-7300, 7300-7600\}^{15}$ Å. The frames were summed within these bands and the raw squared visibility $V_{RAW}$ was derived from the sum. The spectral bands have to be narrow because of the limited coherence of the waves due to the atmospheric turbulence. There are no strong stellar lines in any of the four spectral bands used; the spectral band 7300-7600 Å is affected by the telluric water vapour lines, but even those are smeared out by the low resolution of the spectra.

A calibrator was observed before and after each observation of ξ Tau. Calibrators were selected with the tool Search-Cal (Bonneau et al. 2006), and their list along with their basic properties is given in Table C.5. The instrumental visibility was estimated according to the formulae

$$V_{\text{SCI}}^{2}(u, v) = V_{\text{SCI-RAW}}^{2} \frac{V_{\text{CAL-UD}}^{2}}{V_{\text{CAL-RAW}}^{2}}(u, v), \qquad (\text{C.1})$$

where $V_{\text{SCI}}^2$ is the calibrated visibility of ξ Tau, $V_{\text{SCI-RAW}}^2$ the raw visibility of ξ Tau, $V_{\text{CAL-UD}}^2$ the visibility of a uniform disk with a diameter listed in Table C.5, and $V_{\text{CAL-RAW}}^2$ the raw visibility of a calibrator. To avoid inaccurate observations, we removed all blocks with a $S/N < 2$ and whose optical path delay (OPD)

<footnote n="15">¹⁵ The only difference between the reduction procedure of the observations acquired in 2011 and 2012 is in the choice of the spectral bands. The following bands were used in 2011 $Δλ_{IF(OLD)} = \{5350-5450, 5450-5600\}$ Å, and $Δλ_{IF(OLD)} = \{7000-7200, 7100-7300, 7200-7400\}$ Å.</footnote>

<footnote n="13">¹³ http://gnudatalanguage.sourceforge.net</footnote>

<footnote n="14">¹⁴ http://www.eso.org/~chummel/oyster</footnote>
differed from the mean OPD by more than $2\sigma$. Such blocks usually represent only random noise rather than a physical signal. In rare cases, when the raw visibility of $\xi$ Tau was close to zero, but safely detected, and there was no suitable observation of a calibrator, the raw visibilities of $\xi$ Tau were used in the analysis as if they were calibrated, but they were assigned an error $\Delta V^2 = 0.05$. This allowed us to save more usable observations for very long baselines giving strong constraints by low-visibility measurements.

### C.4. VLTI/AMBER observations

$\xi$ Tau was observed by VLTI/AMBER in 2012 Dec. 03. The observations were acquired in three-telescope mode in J, H, and K bands and the low-resolution regime ($R = 35$). The baselines ranged from 41 m to 139 m.

The unprocessed observations were downloaded from the ESO archive and the reduction was made with the AMBER data reduction software `amdlib`<fnref n="16" />. Following the manual step by step, we applied the bad pixel and flat-field maps, computed pixel-to-visibility matrix, subtracted the dark frame, and performed frame selection based on the fringe S/N ratio (the best 20% were kept).

Four stars were used as calibrators; they are listed along with their basic properties is in Table C.5. The uniform-disk diameters were taken from the JMMC Stellar Diameters Catalogue (Lafrasse et al. 2010). The calibration was made using the library `amdlib`.

Additionally, we had to filter the data. Observations with an effective wavelength at the edges of the J, H, and K bands were not removed, although they had large uncertainties and displayed an abrupt drop in visibility inconsistent with the remaining data. Only observations whose effective wavelength lay in any of following bands $\Delta\lambda \in \{1.155-1.34; 1.49-1.77; 2.02-2.05; 2.075-2.41\}$ µm were used. Furthermore, several observations suddenly had very low visibility compared to neighbouring data, which was very likely caused by an instrumental and/or atmospheric effect. These are data taken at RJD = 56 264.767145, all data with $B/\lambda < 1.76 \times 10^7$, and data taken from RJD = 56 264.776653 to RJD = 56 264.778568 with $B/\lambda > 9.25 \times 10^7$.

### C.5. Night-by-night analysis of NPOI observations

The calibrated visibility and phase estimates are rich enough to permit night-by-night estimation of positions of individual components. Owing to the lower resolution, the NPOI interferometer is almost insensitive to orbit 1 and diameters of the three components (Aa, Ab and B) of $\xi$ Tau. Therefore the system was represented by two point sources, and the relative position of component B and the eclipsing binary was estimated. The results of the night-by-night analysis are given in Table C.1.

The uncertainty ellipses of position of the photocentre of orbit 1 (which is almost identical with its centre of mass because the two components of the eclipsing binary are similar) relative to component B were computed from fits to contours of the $\chi^2$ surfaces at the minima instead of deriving them from the interferometric PSF to take the limitations of fitting a component position very far from the phase centre into account. For the contour we chose 25% above minimum to obtain a reduced $\chi^2$.

<table><caption>**Table C.1.** Astrometric positions of the photocentre of orbit 1 relative to component B derived from night-by-night analysis of MARK III and NPOI observations.</caption><thead><tr><th>RJD</th><th>ρ<br>(mas)</th><th>θ<br>(deg)</th><th>a<br>(mas)</th><th>b<br>(mas)</th><th>α<br>(deg)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-weight: bold;">NPOI</td></tr><tr><td>51 093.906</td><td>9.71</td><td>145.08</td><td>0.831</td><td>0.157</td><td>177.1</td></tr><tr><td>51 097.971</td><td>11.87</td><td>148.00</td><td>0.838</td><td>0.169</td><td>169.1</td></tr><tr><td>51 171.722</td><td>18.31</td><td>327.53</td><td>0.844</td><td>0.153</td><td>173.7</td></tr><tr><td>51 815.933</td><td>7.36</td><td>142.65</td><td>0.842</td><td>0.156</td><td>176.2</td></tr><tr><td>51 835.927</td><td>11.98</td><td>149.63</td><td>0.853</td><td>0.163</td><td>171.7</td></tr><tr><td>52 913.988</td><td>18.72</td><td>327.60</td><td>0.628</td><td>0.111</td><td>151.8</td></tr><tr><td>52 927.944</td><td>18.46</td><td>328.68</td><td>1.962</td><td>0.223</td><td>149.7</td></tr><tr><td>52 930.924</td><td>18.30</td><td>329.24</td><td>0.608</td><td>0.263</td><td>167.0</td></tr><tr><td>55 463.974</td><td>12.48</td><td>148.89</td><td>1.874</td><td>0.256</td><td>152.7</td></tr><tr><td>55 464.970</td><td>12.23</td><td>148.82</td><td>0.675</td><td>0.256</td><td>162.2</td></tr><tr><td>55 465.970</td><td>12.22</td><td>149.96</td><td>0.666</td><td>0.252</td><td>162.2</td></tr><tr><td>55 466.962</td><td>11.74</td><td>149.30</td><td>0.653</td><td>0.254</td><td>162.7</td></tr><tr><td>55 467.963</td><td>11.41</td><td>150.01</td><td>0.651</td><td>0.256</td><td>162.8</td></tr><tr><td>55 468.959</td><td>11.12</td><td>150.03</td><td>0.650</td><td>0.257</td><td>162.7</td></tr><tr><td>55 469.886</td><td>10.93</td><td>150.16</td><td>0.624</td><td>0.274</td><td>180.0</td></tr><tr><td>55 470.955</td><td>10.11</td><td>150.39</td><td>0.643</td><td>0.272</td><td>163.5</td></tr><tr><td>55 999.608</td><td>10.00</td><td>334.30</td><td>2.952</td><td>0.229</td><td>130.3</td></tr><tr><td>56 001.610</td><td>8.30</td><td>335.23</td><td>3.155</td><td>0.250</td><td>126.8</td></tr><tr><td>56 221.917</td><td>5.90</td><td>318.92</td><td>0.424</td><td>0.091</td><td>158.9</td></tr><tr><td>56 227.894</td><td>9.59</td><td>322.81</td><td>0.544</td><td>0.081</td><td>160.8</td></tr><tr><td>56 228.900</td><td>10.52</td><td>324.06</td><td>0.609</td><td>0.098</td><td>154.1</td></tr><tr><td>56 229.901</td><td>11.28</td><td>324.86</td><td>0.620</td><td>0.095</td><td>154.8</td></tr><tr><td>56 230.899</td><td>11.53</td><td>324.83</td><td>0.631</td><td>0.088</td><td>156.9</td></tr><tr><td>56 235.880</td><td>14.12</td><td>325.89</td><td>0.527</td><td>0.081</td><td>160.3</td></tr><tr><td>56 236.878</td><td>14.59</td><td>325.94</td><td>0.497</td><td>0.088</td><td>158.4</td></tr><tr><td>56 237.869</td><td>15.02</td><td>326.35</td><td>0.552</td><td>0.080</td><td>161.1</td></tr><tr><td>56 238.864</td><td>15.45</td><td>326.38</td><td>0.550</td><td>0.080</td><td>161.2</td></tr><tr><td>56 297.679</td><td>4.12</td><td>337.15</td><td>0.787</td><td>0.107</td><td>178.3</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">MARK III</td></tr><tr><td>48 275.689</td><td>18.18</td><td>328.84</td><td>0.852</td><td>0.146</td><td>80.6</td></tr><tr><td>48 548.925</td><td>15.20</td><td>323.62</td><td>2.490</td><td>0.219</td><td>72.8</td></tr><tr><td>48 562.870</td><td>18.02</td><td>327.93</td><td>1.066</td><td>0.314</td><td>85.6</td></tr></tbody></table>

**Notes.** $\rho$ is the separation, $\theta$ the position angle measured counter-clockwise from the north, $a, b$ and $\alpha$ are the semi-major axis, semi-minor axis, and the position angle (again measured counter-clockwise from the north) of the uncertainty ellipse.

An astrometric fit to positions listed in Table C.1 confirms that NPOI is insensitive to the eclipsing binary because the derived orbital properties do not differ significantly from those obtained from a global fit to $V^2$ presented in Table 11.

### C.6. Supplementary materials to the analysis of the spectro-interferometric observations

The spectro-interferometric supplementary material consists of the following tables and figures:

i) Table C.2 lists all available spectro-interferometric observations acquired with the NPOI and MARK III instruments. For each observation the observing stations, its baselines [$B_{\text{min}}; B_{\text{max}}$], phase coverage of orbits 1, and $2\phi_1$, and, $\phi_2$, and associated calibrators are given. The numbering of calibrators is given in Table C.4.

ii) Table C.3 lists all available spectro-interferometric observations acquired with the CHARA/VEGA instrument. For each observation the lengths of the projected baselines $B$

<footnote n="16">¹⁶ The user's manual for the tool `amdlib` was developed at the Jean-Marie Mariotti Center and is publicly available at http://www.jmmc.fr/doc/index.php?search=JMMC-MAN-2720-0001</footnote>
and their orientation $\theta$, the phase coverage of orbits 1, and $2 \phi_1$, and, $\phi_2$, and associated calibrators are given.

iii) Table C.4 lists all calibrators which were used to calibrate the NPOI observations of $\xi$ Tau. For each calibrator its Johnson V magnitude, spectral type, colour index $V-K$, interstellar reddening $E(B-V)$, the minimum amplitude squared visibility $V^2$, and the uniform disk diameter $\theta_{V-K}$ for wavelength range from V to K band are given.

iv) Table C.5 lists all calibrators which were used for calibration of the CHARA/VEGA and VLTI/AMBER

observations. For each calibrator the spectral type, effective temperature $T_{eff}$, gravitational acceleration $log g$, Johnson V (K) magnitude V (K), and the uniform disk diameter in these bands $\theta_V$, $\theta_K$ are given.

v) Figures C.1–C.10 show fits of the global model given by Eq. (9), and corresponding parameters are listed in Table 11.

<table><caption>Table C.2. NPOI and MARK III observations of $\xi$ Tau.</caption><thead><tr><th>RJD<sup>1</sup></th><th>Triangles/baselines</th><th>B<sub>min</sub><br>(m)</th><th>B<sub>max</sub><br>(m)</th><th>φ<sub>1</sub></th><th>φ<sub>2</sub></th><th>Calibrators</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">NPOI</td></tr><tr><td>51 093.906054</td><td>AC-AE-AW</td><td>15</td><td>37</td><td>0.066</td><td>0.897</td><td>09 11 14 18</td></tr><tr><td>51 097.971703</td><td>AC-AE-AW</td><td>16</td><td>37</td><td>0.635</td><td>0.925</td><td>09 11 14 18</td></tr><tr><td>51 171.722612</td><td>AC-AE-AW</td><td>17</td><td>37</td><td>0.954</td><td>0.433</td><td>09 11</td></tr><tr><td>51 815.933549</td><td>AC-AE-AW</td><td>16</td><td>37</td><td>0.096</td><td>0.871</td><td>01 03 09 08</td></tr><tr><td>51 835.927180</td><td>AC-AE-AW</td><td>17</td><td>37</td><td>0.894</td><td>0.008</td><td>20 09 12 14 15 08</td></tr><tr><td>52 913.988131</td><td>AC-AE-W7, AC-AW-W7</td><td>17</td><td>63</td><td>0.742</td><td>0.435</td><td>19</td></tr><tr><td>52 927.944213</td><td>AE-AC, AW-AC</td><td>17</td><td>49</td><td>0.695</td><td>0.531</td><td>05 19</td></tr><tr><td>52 930.924441</td><td>AE-AC, AW-AC</td><td>19</td><td>56</td><td>0.112</td><td>0.551</td><td>05 19</td></tr><tr><td>55 463.974404</td><td>AC-AE, AC-AW</td><td>18</td><td>22</td><td>0.551</td><td>0.000</td><td>06 15 16 17 21</td></tr><tr><td>55 464.970628</td><td>AC-AE, AC-AW</td><td>14</td><td>22</td><td>0.690</td><td>0.007</td><td>06 15 16 17 21</td></tr><tr><td>55 465.970834</td><td>AC-AE, AC-AW</td><td>17</td><td>22</td><td>0.830</td><td>0.014</td><td>06 15 16 17 21</td></tr><tr><td>55 466.962039</td><td>AC-AE, AC-AW</td><td>16</td><td>22</td><td>0.969</td><td>0.021</td><td>06 15 16 17 21</td></tr><tr><td>55 467.963919</td><td>AC-AE, AC-AW</td><td>16</td><td>22</td><td>0.109</td><td>0.028</td><td>06 15 16 17 21</td></tr><tr><td>55 468.959529</td><td>AC-AE, AC-AW</td><td>15</td><td>22</td><td>0.248</td><td>0.034</td><td>06 15 16 17 21</td></tr><tr><td>55 469.886574</td><td>AC-AE, AC-AW</td><td>15</td><td>22</td><td>0.378</td><td>0.041</td><td>06 15 16 17 21</td></tr><tr><td>55 470.955776</td><td>AC-AE, AC-AW</td><td>14</td><td>22</td><td>0.528</td><td>0.048</td><td>06 15 16 17 21</td></tr><tr><td>55 999.608038</td><td>AE-AC, AW-AC</td><td>13</td><td>16</td><td>0.500</td><td>0.690</td><td>10 11 13</td></tr><tr><td>56 001.610251</td><td>AE-AC, AW-AC</td><td>13</td><td>15</td><td>0.780</td><td>0.704</td><td>10 11 13</td></tr><tr><td>56 221.917782</td><td>AC-E6-W7, AC-AE</td><td>19</td><td>73</td><td>0.607</td><td>0.221</td><td>01 07 11</td></tr><tr><td>56 227.894044</td><td>AC-AE-W7, AC-E6-W7</td><td>19</td><td>79</td><td>0.443</td><td>0.262</td><td>07 11 16 02 21</td></tr><tr><td>56 228.900174</td><td>AC-AE-W7</td><td>18</td><td>64</td><td>0.584</td><td>0.269</td><td>07 11 16 02 21</td></tr><tr><td>56 229.901258</td><td>AC-AE-W7</td><td>19</td><td>67</td><td>0.724</td><td>0.276</td><td>07 11 16 02 21</td></tr><tr><td>56 230.899631</td><td>AC-AE-W7, W7-E6</td><td>18</td><td>73</td><td>0.863</td><td>0.283</td><td>07 11 16 02 21</td></tr><tr><td>56 235.880496</td><td>AC-AE-W7, AC-E6-W7</td><td>18</td><td>77</td><td>0.560</td><td>0.317</td><td>07 11 16 02 21</td></tr><tr><td>56 236.878894</td><td>AC-AE-W7, AC-E6-W7</td><td>19</td><td>75</td><td>0.700</td><td>0.324</td><td>07 11 16 02 21</td></tr><tr><td>56 237.869383</td><td>AC-AE-W7, AC-E6-W7</td><td>19</td><td>77</td><td>0.839</td><td>0.331</td><td>07 11 16 02 21</td></tr><tr><td>56 238.864654</td><td>AC-AE-W7, AC-E6-W7</td><td>18</td><td>78</td><td>0.978</td><td>0.338</td><td>07 11 16 04 21</td></tr><tr><td>56 297.679445</td><td>AC-AE-W6</td><td>0</td><td>53</td><td>0.208</td><td>0.743</td><td>10 11</td></tr><tr><td colspan="7" style="text-align: center; font-weight: bold;">Mark III</td></tr><tr><td>48 275.689</td><td>NF-SF</td><td>29</td><td>30</td><td>0.725</td><td>0.484</td><td>μ Tau, η Tau</td></tr><tr><td>48 548.925</td><td>ND-SC</td><td>14</td><td>15</td><td>0.958</td><td>0.366</td><td>μ Tau, η Tau</td></tr><tr><td>48 562.870</td><td>ND-SC</td><td>14</td><td>15</td><td>0.909</td><td>0.462</td><td>μ Tau, η Tau</td></tr></tbody></table>

**Notes.** Column 2 lists the configuration used as a triple of stations (e.g. “AC-AE-W7”, using astrometric stations centre and east, as well as imaging station W7) if data from all three involved baselines were recorded, including the corresponding closure phase. When a single baseline is listed, only squared visibility data were recorded but no closure phases on that baseline. $\phi_1$ ($\phi_2$) refers to mean orbital phase of orbit 1 (2) at the given RJD. Calibrator numbers correspond to the numbering in Table C.4. (1) Mean RJD. Ephemeris used to compute the orbital phases: $1 - T_{min,1}(RJD) = 7.1467 \times E + 56 224.7246, 2 - T_{p,2}(RJD) = 145.17 \times E + 55 609.05$, where E is the epoch, $T_{min,1}$ the epoch of the primary minimum of orbit 1, $T_{p,2}$ the epoch of the periastron passage of orbit 2.
<table><caption>**Table C.3.** Journal of CHARA/VEGA spectro-interferometric observations of ξ Tau.</caption><thead><tr><th>RJD</th><th>NB</th><th>φ<sub>1</sub></th><th>φ<sub>2</sub></th><th>B (m)</th><th>θ (deg)</th><th>Cal.</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center;">The 2011 run</td></tr><tr><td rowspan="3">55 825.8907</td><td>3-1</td><td rowspan="3">0.193</td><td rowspan="3">0.488</td><td>064.6</td><td>-155.9</td><td rowspan="3">1, 2</td></tr><tr><td>3-2</td><td>150.0</td><td>-160.6</td></tr><tr><td>3-3</td><td>217.2</td><td>-158.9</td></tr><tr><td rowspan="3">55 846.8703</td><td>3-1</td><td rowspan="3">0.129</td><td rowspan="3">0.633</td><td>065.8</td><td>-154.2</td><td rowspan="3">2</td></tr><tr><td>3-2</td><td>155.8</td><td>-159.3</td></tr><tr><td>3-3</td><td>221.4</td><td>-157.8</td></tr><tr><td rowspan="3">55 850.8130</td><td>3-1</td><td rowspan="3">0.680</td><td rowspan="3">0.660</td><td>147.0</td><td>-160.9</td><td rowspan="3">1, 2</td></tr><tr><td>3-2</td><td>154.1</td><td>-090.1</td></tr><tr><td>3-3</td><td>244.8</td><td>+056.0</td></tr><tr><td rowspan="3">55 854.8645</td><td>3-1</td><td rowspan="3">0.247</td><td rowspan="3">0.688</td><td>065.6</td><td>-153.1</td><td rowspan="3">1, 2</td></tr><tr><td>3-2</td><td>156.2</td><td>-158.3</td></tr><tr><td>3-3</td><td>221.6</td><td>-156.8</td></tr><tr><td rowspan="3">55 854.9548</td><td>3-1</td><td rowspan="3">0.260</td><td rowspan="3">0.688</td><td>135.3</td><td>-148.5</td><td rowspan="3">1, 2</td></tr><tr><td>3-2</td><td>172.7</td><td>-057.7</td></tr><tr><td>3-3</td><td>217.7</td><td>+084.1</td></tr><tr><td rowspan="3">55 856.8928</td><td>3-1</td><td rowspan="3">0.531</td><td rowspan="3">0.702</td><td>063.3</td><td>-149.9</td><td rowspan="3">2, 3</td></tr><tr><td>3-2</td><td>152.5</td><td>-155.6</td></tr><tr><td>3-3</td><td>216.3</td><td>-154.2</td></tr><tr><td colspan="7" style="text-align: center;">The 2012 run</td></tr><tr><td>56 194.8118</td><td>2-1</td><td>0.814</td><td>0.029</td><td>054.8</td><td>-156.0</td><td>1</td></tr><tr><td>56 194.9180</td><td>2-1</td><td>0.829</td><td>0.030</td><td>065.9</td><td>-154.1</td><td>1</td></tr><tr><td rowspan="3">56 197.8894</td><td>3-1</td><td rowspan="3">0.245</td><td rowspan="3">0.050</td><td>065.6</td><td>-155.2</td><td rowspan="3">1</td></tr><tr><td>3-2</td><td>153.5</td><td>-160.1</td></tr><tr><td>3-3</td><td>218.7</td><td>-158.6</td></tr><tr><td rowspan="3">56 197.9362</td><td>3-1</td><td rowspan="3">0.252</td><td rowspan="3">0.051</td><td>065.0</td><td>-152.1</td><td rowspan="3">1</td></tr><tr><td>3-2</td><td>155.6</td><td>-157.5</td></tr><tr><td>3-3</td><td>220.4</td><td>-155.9</td></tr><tr><td rowspan="3">56 200.0052</td><td>3-1</td><td rowspan="3">0.541</td><td rowspan="3">0.065</td><td>106.3</td><td>-002.5</td><td rowspan="3">1</td></tr><tr><td>3-2</td><td>203.9</td><td>-060.1</td></tr><tr><td>3-3</td><td>276.1</td><td>-041.1</td></tr><tr><td rowspan="3">56 200.0306</td><td>3-1</td><td rowspan="3">0.545</td><td rowspan="3">0.065</td><td>099.2</td><td>-000.3</td><td rowspan="3">1</td></tr><tr><td>3-2</td><td>207.7</td><td>-056.8</td></tr><tr><td>3-3</td><td>278.3</td><td>-039.2</td></tr><tr><td>56 226.9927</td><td>2-1</td><td>0.317</td><td>0.251</td><td>045.4</td><td>-125.8</td><td>1</td></tr><tr><td>56 227.0299</td><td>2-1</td><td>0.323</td><td>0.251</td><td>040.0</td><td>-109.9</td><td>1</td></tr><tr><td>56 227.8758</td><td>2-1</td><td>0.441</td><td>0.257</td><td>031.3</td><td>+100.3</td><td>1</td></tr><tr><td>56 227.9720</td><td>2-1</td><td>0.454</td><td>0.258</td><td>033.4</td><td>+117.0</td><td>1</td></tr></tbody></table>

**Notes.** $\phi_1$ ($\phi_2$) denotes the orbital phase of orbit 1 (2), B the mean length of the projected baseline, $\theta$ the position angle of the projected baseline. The calibrator stars are identified as follows: 1 – HD 21686, 2 – HD 18604, and 3 – HD 26793. Ephemeris used to compute the orbital phases: $1 - T_{\text{min},1}(\text{RJD}) = 7.1467 \times E + 56.224.7246$, $2 - T_{\text{p},2}(\text{RJD}) = 145.17 \times E + 55.609.05$, where E is the epoch, $T_{\text{min}}^j$ the epoch of the primary minimum of orbit 1, $T_{\text{p},2}$ the epoch of the periastron passage of orbit 2.

<table><caption>**Table C.4.** List of NPOI calibrators used for ξ Tau.</caption><thead><tr><th>No</th><th>HD</th><th>Type</th><th>V</th><th>V-K</th><th>E(B-V)</th><th>V<sup>2</sup></th><th>θ<sub>V-K</sub> (mas)</th></tr></thead><tbody><tr><td>01</td><td>886</td><td>B2IV</td><td>2.83</td><td>-0.940</td><td>0.010</td><td>0.85</td><td>0.498</td></tr><tr><td>02</td><td>7804</td><td>A3V</td><td>5.16</td><td>0.239</td><td>-0.010</td><td>0.91</td><td>0.366</td></tr><tr><td>03</td><td>7964</td><td>A3V</td><td>4.76</td><td>0.224</td><td>-0.050</td><td>0.97</td><td>0.434</td></tr><tr><td>04</td><td>11171</td><td>F3III</td><td>4.65</td><td>0.778</td><td>-0.035</td><td>0.73</td><td>0.653</td></tr><tr><td>05</td><td>12216</td><td>A2V</td><td>3.98</td><td>0.059</td><td>-0.060</td><td>0.98</td><td>0.562</td></tr><tr><td>06</td><td>16582</td><td>B2IV</td><td>4.07</td><td>-0.632</td><td>0.020</td><td>0.99</td><td>0.343</td></tr><tr><td>07</td><td>17081</td><td>B7V</td><td>4.25</td><td>-0.255</td><td>-0.010</td><td>0.89</td><td>0.403</td></tr><tr><td>08</td><td>20630</td><td>G5Vvar</td><td>4.83</td><td>1.873</td><td>0.000</td><td>0.86</td><td>0.956</td></tr><tr><td>09</td><td>23630</td><td>B7III</td><td>2.90</td><td>0.264</td><td>0.010</td><td>0.85</td><td>0.981</td></tr><tr><td>10</td><td>24760</td><td>B0.5V</td><td>2.88</td><td>-0.833</td><td>0.110</td><td>0.91</td><td>0.519</td></tr><tr><td>11</td><td>25490</td><td>A1V</td><td>3.91</td><td>0.127</td><td>0.020</td><td>0.77</td><td>0.600</td></tr><tr><td>12</td><td>37128</td><td>B0Ia</td><td>1.70</td><td>-0.573</td><td>0.040</td><td>0.84</td><td>1.012</td></tr><tr><td>13</td><td>76756</td><td>A5m</td><td>4.20</td><td>0.256</td><td>0.190</td><td>0.98</td><td>0.582</td></tr><tr><td>14</td><td>184006</td><td>A5Vn</td><td>3.79</td><td>0.192</td><td>-0.010</td><td>0.93</td><td>0.668</td></tr><tr><td>15</td><td>192696</td><td>A3IV-Vn</td><td>4.30</td><td>0.222</td><td>0.030</td><td>0.96</td><td>0.536</td></tr><tr><td>16</td><td>195810</td><td>B6III</td><td>4.03</td><td>-0.351</td><td>0.020</td><td>0.88</td><td>0.421</td></tr><tr><td>17</td><td>196724</td><td>A0V</td><td>4.82</td><td>-0.034</td><td>0.000</td><td>0.99</td><td>0.360</td></tr><tr><td>18</td><td>213558</td><td>A1V</td><td>3.77</td><td>-0.081</td><td>0.000</td><td>0.95</td><td>0.568</td></tr><tr><td>19</td><td>214923</td><td>B8.5V</td><td>3.40</td><td>-0.166</td><td>0.003</td><td>0.85</td><td>0.635</td></tr><tr><td>20</td><td>216735</td><td>A1V</td><td>4.90</td><td>0.060</td><td>-0.010</td><td>0.99</td><td>0.366</td></tr><tr><td>21</td><td>217891</td><td>B6Ve</td><td>4.53</td><td>-0.220</td><td>0.030</td><td>0.92</td><td>0.360</td></tr></tbody></table>

**Notes.** V(K) is the apparent magnitude in the Johnson V(K) filter, E(B-V) the interstellar reddening, V$^2$ the squared visibility, and $\theta$ the uniform disk diameter.
<table><caption>**Table C.5.** List of stars used for calibration of CHARA/VEGA and VLTI/AMBER observations.</caption><thead><tr><th rowspan="2">Calibrator</th><th rowspan="2"><i>T</i><sub>eff</sub><br>(K)</th><th rowspan="2">log <i>g</i>_{[cgs]}</th><th rowspan="2"><i>V</i><br>(mag)</th><th rowspan="2"><i>K</i><br>(mag)</th><th colspan="3">Parameter</th><th rowspan="2">θ<sub>H</sub><br>(mas)</th><th rowspan="2">θ<sub>K</sub><br>(mas)</th></tr><tr><th>θ<sub>V</sub><br>(mas)</th><th>θ<sub>R</sub><br>(mas)</th><th>θ<sub>J</sub><br>(mas)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">CHARA/VEGA</td></tr><tr><td>HD 21686</td><td>9790</td><td>4.1</td><td>5.125</td><td>5.167</td><td>0.245(18)</td><td>0.247(18)</td><td>0.251(18)</td><td>0.252(18)</td><td>0.252(18)</td></tr><tr><td>HD 18604</td><td>13 000</td><td>3.4</td><td>4.703</td><td>4.910</td><td>0.257(18)</td><td>0.257(18)</td><td>0.262(18)</td><td>0.262(18)</td><td>0.263(18)</td></tr><tr><td>HD 26793</td><td>10 500</td><td>4.0</td><td>5.210</td><td>5.357</td><td>0.207(15)</td><td>0.209(15)</td><td>0.212(15)</td><td>0.213(15)</td><td>0.214(15)</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">VLTI/AMBER</td></tr><tr><td>HD 25490</td><td>9500</td><td>4.1</td><td>3.891</td><td>3.783</td><td>0.513(37)</td><td>0.518(37)</td><td>0.526(37)</td><td>0.529(37)</td><td>0.530(37)</td></tr><tr><td>HD 34909</td><td>4660</td><td>2.1</td><td>7.987</td><td>5.775</td><td>0.310(23)</td><td>0.315(23)</td><td>0.323(23)</td><td>0.326(23)</td><td>0.328(23)</td></tr><tr><td>HD 38277</td><td>5700</td><td>4.4</td><td>7.119</td><td>5.597</td><td>0.318(23)</td><td>0.322(23)</td><td>0.329(23)</td><td>0.331(23)</td><td>0.333(23)</td></tr><tr><td>HD 38406</td><td>5790</td><td>4.4</td><td>8.197</td><td>6.735</td><td>0.186(14)</td><td>0.188(14)</td><td>0.192(14)</td><td>0.193(14)</td><td>0.194(14)</td></tr></tbody></table>

**Notes.** All data were taken from the JMMC Stellar Diameters Catalogue (Lafrasse et al. 2010). *T*<sub>eff</sub> denotes the effective temperature, *g* the gravitational acceleration, *θ*<sub>X</sub> the uniform-disk diameter in the passband X, *X* the magnitude in the passband X.
<figure><img src="image_19.png" /><figcaption>**Fig. C.1.** Best-fitting model (part one) plotted against the observations from the MARKIII and NPOI spectro-interferometers. In each panel, the observed squared visibility $V^2$ is plotted with red triangles; the model corresponding to parameters listed in Table 11 is denoted with black points. Residuals of the fit are shown below each panel. The mean acquisition date, the corresponding mean heliocentric Julian date, and the instrument are indicated above each panel.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. C.2.** Best-fitting model (part two) plotted against the observations from the NPOI spectro-interferometer. In each panel, the observed squared visibility $V^2$ is plotted with red triangles; the model corresponding to parameters listed in Table 11 is denoted with black points. Residuals of the fit are shown below each panel. The mean acquisition date, the corresponding mean reduced heliocentric Julian date, and the instrument are indicated above each panel.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

**Fig. C.3.** Best-fitting model (part three) plotted against the observations from the NPOI spectro-interferometer. In each panel, the observed squared visibility $V^2$ is plotted with red triangles; the model corresponding to parameters listed in Table 11 is denoted with black points. Residuals of the fit are shown below each panel. The mean acquisition date, the corresponding mean reduced heliocentric Julian date, and the instrument are indicated above each panel.
<figure><img src="image_22.png" /><figcaption>**Fig. C.4.** Best-fitting model (part four) plotted against the observations from the NPOI spectro-interferometer. In each panel, the observed closure phase $T_3\phi$ is plotted with red triangles; the model corresponding to parameters listed in Table 11 is denoted with black points. Residuals of the fit are shown below each panel. The mean acquisition date, the corresponding mean reduced heliocentric Julian date, and the instrument are indicated above each panel.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. C.5.** Best-fitting model (part five) plotted against the observations from the NPOI spectro-interferometer. In each panel, the observed closure phase $T_3\phi$ is plotted with red triangles; the model corresponding to parameters listed in Table 11 is denoted with black points. Residuals of the fit are shown below each panel. The mean acquisition date, the corresponding mean reduced heliocentric Julian date, and the instrument are indicated above each panel.</figcaption></figure>