© The Authors 2026

1. Introduction

Triple, quadruple, or quintuple stellar systems are not just ‘more complicated binaries’ (Tokovinin & Moe 2020; Tokovinin 2021; Borkovits et al. 2020; Borkovits 2022). They were not assembled randomly, but their varying architectures and occurrences (2+1, 3+1, 2+2, 3+2, etc.) suggest several evolutionary processes, for example, migration, Kozai cycles, tidal friction, magnetic friction, envelope expansion, mass transfer, close encounters, or ejections.

Less complicated binaries form a foundation of the universal distance scale (Pietrzyński et al. 2013, 2019; Gallenne et al. 2023), provided their properties were measured with percentage-level precision. Multiples allowed us to achieve even higher precision due to additional dynamical constraints (Borkovits 2022; Brož et al. 2025). Two or more orbits might also serve as an additional verification.

The ξ Tau system is composed of four components (Aa, Ab, B, C). From a Solar-System perspective, the inner eclipsing binary has a semi-major axis of only 0.11 au, the third component is similar to the Earth (∼1.1 au), and the fourth component is similar to Neptune (29 au); the total mass is over 9 M_⊙.

In our previous work (Nemravová et al. 2016), we used space-based photometry from the Microvariability and Oscillations of Stars Satellite (MOST) (Walker et al. 2003), which allowed a clear detection of mutual perturbations; for example, the eclipse-timing variations (ETVs) exhibit an amplitude of about half an hour. In this work, we aim to use space-based photometry from the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al. 2015) and an additional set of MOST observations from 2017, supplemented by ground-based spectroscopy from Cerro Tololo Inter-American Observatory (CTIO) (Tokovinin et al. 2013). Consequently, the high-accuracy observations extend over a substantially longer time span. Their analysis should allow a much better characterization of the system parameters, including evidence of the evolution of ξ Tau orbital architecture on various timescales (from days to decades).

A related scientific question is whether or not this stellar system also contains exoplanets. More than 800 exoplanets are already known in binaries (Thebault & Haghighipour 2015; Thebault & Bonanni 2025). More than 30 are circumbinary, where a secondary is commonly located at 0.1 au and an exoplanet at ∼0.5 au, which is close to the stability limit (Holman & Wiegert 1999). This rough statistics is biased, though, especially for O or B stars, which are much brighter than exoplanets. One should keep this in mind when studying bright stellar systems.

Hereinafter, we describe only new observations of ξ Tau (Sect. 2), since old observations were described in Nemravová et al. (2016). Subsequently, we describe models of ξ Tau of various complexities, the reference one (Sect. 3.1), the interferometric one (Sect. 3.2), the all-data model (Sect. 3.3), the model with tides (Sect. 3.4), and the model with five components (Sect. 3.5). The temporal evolution of the ξ Tau system is discussed in Sect. 4 and additional aspects (parallax, oscillations) are detailed in Sect. 5. The conclusions are presented in Sect. 6.

2. New observations

2.1. Hvar photometry

Systematic photometry was collected at the Hvar Observatory with the 0.65-m reflector and photoelectric photometer, initially in UBV bands (2007–2013) and later in UBVR bands (2013–2025). After a recent revision of 52 years of Hvar observations (Božić et al. 2026), all-sky magnitudes of all comparison stars were derived relatively to the Johnson standards and all observations were carefully transformed to the standard Johnson UBV system (Harmanec et al. 1994). The latest reduction includes also temporal variation of the extinction coefficients in the course of night. The Hvar mean values for comparison stars, 4 Tau = HD 21686, and 6 Tau = HD 21933, which were added to the differential magnitudes, are V = 5.150, B − V = −0.045, U − B = −0.093, V − R = −0.016; V = 5.774, B − V = −0.076, U − B = −0.300, V − R = −0.025.

We note that the transparency curve of the R filter closely corresponds to that of the standard Cousins R filter. However, because we were not able to find enough northern bright standard stars with the Cousins R_c values, we derived robust mean values of Johnson V − R indices from Johnson et al. (1966) and reduced our observations to the Johnson R magnitude.

2.2. MOST photometry

Broadband photometry obtained by the MOST satellite (Walker et al. 2003) consists of two separate sets of observations. Back in 2012, J. Nemravová and P. Harmanec submitted a successful application for observations of ξ Tau to the MOST allocation committee, and continuous light curve observations were indeed obtained. These observations were published in Nemravová et al. (2016), along with the discovery of low-amplitude photometric oscillations on the time scale of about 0.4 d. We use these data also here.

The second data set was obtained on a commercial basis at the end of 2017 and represents the last set of scientifically usable MOST data before the communication with the satellite was lost forever. The data were downloaded at the Vienna tracking station and their initial reduction was carried out by T. Kallinger.

2.3. TESS photometry

Photometry obtained by the TESS satellite (Ricker et al. 2015) was also included. TESS records red optical light, with a wide bandpass spanning roughly 600−1000 nm centered on the traditional Cousins I band. TESS observed ξ Tau in six sectors (31, 42, 43, 44, 70, 71). Basic reductions and cleaning were performed by Jon Labadie-Bartz. Light curves were extracted from the full frame images with the Lightkurve package (Cardoso et al. 2018) using simple aperture photometry, including the saturated columns and their ‘spillover’ pixel end caps, and the non-saturated halo pixels, as is standard for moderately saturated bright stars. The observing cadence was ten minutes (sectors 31, 42, 43, 44) and 200 seconds (sectors 70, 71). Background subtraction was used as the preferred detrending method. ξ Tau is fairly isolated on the sky and so blending from neighboring stars did not contribute to the observed signals.

The series of light curves represents the key observational constraint, because the exact epochs of eclipses are sensitive to the perturbations by the component B. These include both the geometric effects due to finite speed of light and variations in position of the eclipsed star with respect to the observer (usually termed light-time effect), and variations in the orbital elements of the eclipsing binary due to the gravitational perturbations by the component B (often termed physical timing variations).

2.4. CTIO/CHIRON spectroscopy

The CHIRON instrument (Tokovinin et al. 2013) at the CTIO 1.5-m telescope is an echelle spectrograph. Its resolution is R = 79000 (higher resolution was not used), spectral range 410 to 890 nm, and efficiency at least 6%. It is equipped with the back-illuminated detector CCD231-84 Teledyne E2v, with 4094 × 4112 pixels and gain of η = 1.3 adu/e⁻. We performed a standard reduction of spectra by applying a dark frame and a flat field and by calibrating the wavelength scale with a ThAr lamp spectrum. For a detailed description of rectification, see Appendix A.

In the course of observational campaign, we obtained 277 spectra of ξ Tau (excluding defective ones). The spectral range was limited to 450 to 890 nm, divided into 62 orders (including some overlaps and gaps), with 801 pixels per order. The exposure time was 60 s, the peak signal reaching S = 160000 adu, resulting in the signal-to-noise ratio $S/N=\sqrt{S/\eta }=350$. We observed sequences of five spectra. The total time span was 129 days, from Sep. 15 2021 to Jan. 21 2022, covering the short orbit (Aa+Ab) 18 times and about 90% of the long orbit ((Aa+Ab)+B).

2.5. Gaia parallax

The parallax in the Gaia EDR3/DR3 release (Bailer-Jones et al. 2021; Vallenari et al. 2023) is π = (16.8 ± 0.7) mas, which is equivalent to d = (59.6 ± 2.4) pc. On the contrary, the previous modeling of Nemravová et al. (2016) lead to the different distance d = (67.9 ± 1.0) pc, which is likely the correct one, because ξ Tau exhibits photocentre motion, negatively impacting the parallactic measurement.

2.6. WDS astrometry

Seven new astrometric measurements from the Washington Double Star (WDS) catalogue (Mason et al. 2001; Tokovinin et al. 2020) were included. Together with a removal of two outlier points (in Julian date, TDB) 2456314.597291, 2456936.276038, the new datset provides a slightly different outer orbit solution, compared to the previous modeling in Nemravová et al. (2016).

2.7. Other data

A description of other data was presented in Nemravová et al. (2016), namely of old radial velocities (RVs), eclipse timing variations (ETVs), eclipse durations, interferometric visibility, closure phase, triple product, and the spectral-energy distribution (SED). Of course, these data sets were also used in this study in order to fully constrain the model of ξ Tau.

3. Configuration of ξ Tau

To describe the dynamics of ξ Tau, we used the eponymous model Xitau¹ (Nemravová et al. 2016; Brož 2017; Brož et al. 2021; Marchis et al. 2021; Brož et al. 2022a,b; Ferrais et al. 2022; Brož et al. 2023; Fuksa et al. 2023; Oplištilová et al. 2023). It accounts for mutual, N-body perturbations, relativistic effects f_ppnⁱ in the parametrised post-Newtonian (PPN) approximation (Standish & Williams 2006), oblateness f_oblatⁱ (Kaula 1966), and tides f_tideⁱ (Mignard 1979). The corresponding equation of motion for the position vector, r_i, of the i-th component in the barycentric reference frame (r_ij ≡ r_i − r_j, m_j mass of the j-th component),

$$\begin{array}{c}\hfill {\ddot{\mathit{r}}}_i=-{\displaystyle \sum _{j\ne i}}\frac{G{m}_j}{r_{\mathit{ij}}^3}{\mathit{r}}_{\mathit{ij}}+{\mathit{f}}_{\mathrm{ppn}}^i+{\mathit{f}}_{\mathrm{oblat}}^i+{\mathit{f}}_{\mathrm{tide}}^i,\end{array}$$(1)

is integrated numerically, using the Bulirch-Stoer algorithm (Bulirsch & Stoer 1966; Levison & Duncan 1994), with adaptive time stepping and the relative precision ε = 10⁻⁸. Since the system has four components, we have a total of 48 dynamical and radiative parameters (see Table 1 notes). For the minimisation of the χ², we used the simplex (Nelder & Mead 1965) or subplex (Rowan 1990) algorithms.

To constrain the model, one should use all kinds of observations by minimising the metric

$$\begin{array}{cc}\hfill {\chi }^2=& {\chi }_{\mathrm{sky}}^2+{\chi }_{\mathrm{rv}}^2+{\chi }_{\mathrm{etv}}^2+{\chi }_{\mathrm{ecl}}^2+{\chi }_{\mathrm{vis}}^2+{\chi }_{\mathrm{clo}}^2+{\chi }_{\mathrm{t}3}^2+{\chi }_{\mathrm{lc}}^2\hfill \\ \hfill & +{\chi }_{\mathrm{syn}}^2+{\chi }_{\mathrm{sed}}^2+{\chi }_{\mathrm{sed}2}^2,\hfill \end{array}$$(2)

where individual terms are standard O − Cs squared, divided by uncertainties σ squared, for astrometry (SKY), radial velocities (RV), minima timings (ETV), eclipse durations (ECL), visibility (VIS), closure phase (CLO), triple product (T3), light curves (LC), normalised spectra (SYN), spectral-energy distribution (SED), and relative brightness of the fourth component (SED2).

It is, however, often impossible to converge all parameters at once. We thus proceeded sequentially, by subsequently minimising the following metrics:

$$\begin{array}{cc}& {\chi }^2={\chi }_{\mathrm{etv}}^2+{\chi }_{\mathrm{ecl}}^2,\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\chi }^2={\chi }_{\mathrm{sky}}^2+{\chi }_{\mathrm{rv}}^2+{\chi }_{\mathrm{etv}}^2+{\chi }_{\mathrm{ecl}}^2+{\chi }_{\mathrm{sed}}^2+{\chi }_{\mathrm{sed}2}^2,\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\chi }^2={\chi }_{\mathrm{sky}}^2+{\chi }_{\mathrm{rv}}^2+{\chi }_{\mathrm{etv}}^2+{\chi }_{\mathrm{ecl}}^2+{\chi }_{\mathrm{vis}}^2+{\chi }_{\mathrm{clo}}^2+{\chi }_{\mathrm{t}3}^2+{\chi }_{\mathrm{sed}}^2+{\chi }_{\mathrm{sed}2}^2,\hfill \end{array}$$(5)

always carefully considering which parameters should be free versus fixed. For instance, when using Eq. (3), all radiative parameters must be fixed because they would be totally unconstrained. Also, when some parameters are correlated, such as the total mass, m_sum, and the distance, d (Kepler 1619), it is very useful to optimise grids of models in order to find the global best-fit solution and eodem tempore exclude other solutions.

3.1. Reference model ((Aa+Ab)+B)+C

The reference model was adopted from Nemravová et al. (2016). We added new observations, though, in a simplified form, i.e. RVs instead of full spectra and minima timings instead of full light curves. This allowed us to first focus on the dynamical (not radiative) parameters. We also used a slightly different set of parameters (P instead of a, log e instead of e, ϖ instead of ω, λ instead of M).

3.1.1. Mirror solutions.

We examined eight mirror solutions with ETVs, which are very sensitive to mutual perturbations. We also used eclipse durations to constrain the precession. Specifically, for the inclinations i₁, i₂ and the longitudes of nodes Ω₁, Ω₂, we tested the following combinations of discrete values 86.5, 87.0, 93.0, 93.5° and 147.5, 148.5, 327.5, 328.5°, respectively. Other orbital parameters were free, except for the total mass, m_sum. There is only one combination that explains the ETVs (Table 1). Their overall amplitude with respect to a two-body model is about 0.02 d (Nemravová et al. 2016), although the residuals with respect to the N-body model are < 0.002 d (Fig. 2). The observed eclipse duration decreases from 0.27 to 0.25 d, and the observed eclipse depth decreases from 0.110 to 0.077 mag. However, ETVs are not enough to obtain all dynamical parameters. We thus considered six data sets (SKY, RV, ETV, ECL, SED, SED2), re-converged the reference model, and obtained the free parameter values in Table 1 and the dependent parameter values in Table C.2. The unreduced χ_etv², χ_ecl² values are already comparable to the number of degrees of freedom, ν = N − M; on the other hand, χ_rv² is relatively larger, due to remaining systematics in RVs of the order of less than 1 km s⁻¹. This, however, might be related to systematics of RVs, which will be circumvented by using full spectra.

Fig. 1. Rectification of a CTIO/CHIRON echelle spectrum (ktc00027) in the region of the Hα line (6563 Å), using the resimplex3 method. Top: Synthetic spectra were computed from our model in Table 1, ‘all-data’ column. Four component spectra (Aa, Ab, B, C) and the resulting spectrum are plotted. Second row: Observed spectrum from CHIRON (black), orders r = 39 to 41, and the corresponding blaze function (green). Third row: Rectified spectrum (black), its comparison to the synthetic spectrum (blue), and corresponding residuals (red). The rectification was done with respect to the continuum level, with interpolation of the blaze parameters over wide lines (cyan) and problematic orders, also avoiding locations of the telluric lines (green). Bottom: Residuals plotted separately.

Fig. 2. Reference model of ξ Tau, showing ETVs (top), eclipse durations (middle), perturbations with respect to a two-body ephemeris (bottom), MOST and TESS observations (blue), model (grey), residuals (red), and, for reference, observed minima (green). The total χ2 = 1553 (see Table 1 for individual contributions). The mean ephemeris used in the bottom panel is P = 7.146665 d, T0 = 2456224.704705. Our model is certainly able to explain the major perturbations, which are substantial (±0.02 d, or ±30 min), but some systematics remain to be explained (up to 0.002 d, or 3 min).

3.2. Models constrained by interferometry

Next, we computed a grid of models including three more interferometric data sets (VIS, CLO, T3). Our aim was to constrain the total mass, m_sum, and the distance, d, and to exclude alternative solutions. We used 121 pairs of m_sum and d values and fixed these parameters, but all other parameters were set free; each model was re-converged so that it adapts ‘at all cost’.

The results in Fig. 3 show logical correlations between m_sum, d, as required by astrometric or interferometric data sets (i.e. angular separations of the third and fourth components). However, other data sets (RV, ETV, SED, SED2) are orthogonal to this because some scales are set absolutely, in km s⁻¹, days, and W m⁻² m⁻¹, respectively. This implies a uniqueness of the best-fit solution (Table 1). How this solution with χ² = 1553 looks graphically is shown in Fig. 4. In either case, the ξ Tau distance must be around 66.0 pc; the Gaia distance of 59.6 pc was excluded.

Fig. 3. Best-fit models of ξ Tau constrained by multi-technique observations for two fixed parameters, the total mass msum and the distance d. The χ2 values (colours) indicate best fit (cyan), good fits (blue), and poor fits (orange); the overall best fit (red) has χ2 = 1661 (unreduced). We tested 121 pairs of parameters and 3000 sequential iterations of simplex, i.e. 3.6 × 105 models in total. Other dynamical parameters (except msum, d) were free; other radiative parameters were fixed (Nemravová et al. 2016). We constrained the model by SKY, RV, ETV, ECL, SED, SED2 data sets; we did not use interferometry (but we used SKY), spectroscopy (but RV), or light curves (but ETV, ECL). We assumed unit weights, except wetv = wecl = 100, because we wanted to exclude all models with incorrect dynamics. The two parameters are often positively correlated (e.g. SKY, VIS, CLO, T3), as predicted (Kepler 1619), but some data sets are distance independent or orthogonal (RV, ETV, SED, SED2), and the model is considered well constrained. After re-convergence (green), the χ2 further decreased to 1553. The Gaia distance of 59.6 pc was excluded. Additionally, VIS, CLO, and T3 were also computed, but their weights were set to zero, otherwise our model would exhibit tension between VIS and CLO.

Fig. 4. Best-fit model from Fig. 3 with χ2 = 1553 (unreduced). Plots for each data set, SKY, RV, ETV, ECL, VIS, CLO, T3, SED, and SED2, contain observations (blue), synthetic data (yellow), and residuals (red). The best-fit distance is d = 66.0 pc. One can notice some systematics for ETV (see MOST 2017), VIS, or CLO data sets.

3.2.1. MCMC analysis.

For this model, we estimated the uncertainties of model parameters with the Markov chain Monte Carlo (MCMC) method (Markov 1906), which should reflect the uncertainties of observational data. Here, we used all dynamical parameters, but not all radiative ones, in order to prevent systematics from affecting the resulting uncertainties. If any systematics are present (e.g. tension between VIS and CLO) it will result in very small, unrealistic uncertainties. The standard corner plot is shown in Fig. B.4.

3.3. Models constrained by CHIRON and TESS observations

Next, we included 44 high-resolution spectra (from CHIRON, one per night) and eight high-precision light curves (from MOST and TESS) and computed another grid of models. Our aim was to better constrain the radiative parameters, in particular the effective temperatures T₁ and T₃. We used 121 pairs of T₁ and T₃ values and fixed them, but all other parameters were left free.

The results in Fig. 5 show some data sets that are insensitive to the radiative parameters (e.g. SKY, RV, ETV, ECL), but others constrain the temperatures very well (VIS, LC, SYN). In particular, light curves and eclipse depths are very ‘strict’ constraints (T₁ ≃ 11000 K, T₃ ≃ 14000 K), which is fully in agreement with spectral line profiles. The third component does not eclipse anything, but its brightness contributes as a third light. The line profiles of Aa, Ab, B components are blended, but easy to disentangle because of different rotation speeds (Fig. B.2).

Fig. 5. Same as Fig. 3 but for two different fixed parameters, the effective temperatures T1and T3. This time, light curves (LC) and synthetic spectra (SYN) were included. The overall best fit χ2 = 3639180 (unreduced). One can notice some systematics; for example, there is tension between VIS and CLO, the best ECL is offset, the best SED is offset, and the best SED2 is very offset.

3.3.1. Tension.

Looking at individual contributions to the χ² versus the position of the best-fit model (Fig. 5, red circle), one can notice the following problems. On one hand, the CLO data set exhibits tension (with respect to VIS), where better fits of CLO data require lower T₃. On the other hand, SED, SED2 data sets exhibit different tension (with regard to LC, SYN), requiring higher T₃. It is clearly not easy to correct these issues by changing T₃, or, as a matter of fact, by changing any other parameter; we note that all other parameters were free.

3.3.2. Systematics.

There might be various reasons for such inconsistencies. First, our model could be internally inconsistent; for example, the light curve algorithm (Wilson & Devinney 1971; Wilson et al. 2010), which is technically comparable to Prša et al. (2016), could be incompatible with grids of synthetic spectra (Lanz & Hubený 2003, 2007; Palacios et al. 2010; de Laverny et al. 2012; Husser et al. 2013) used for fitting of the SYN data set. Second, some approximations might be inadequate; for example, the third component is a fast-rotating star, but we used a synthetic spectrum suitable for slow-rotating stars. The same is true for the interferometric quantities, or VIS, CLO, and T3, where we assumed limb-darkened discs (Hanbury Brown et al. 1974). Third, some fixed parameters might be set incorrectly, for example the linear limb-darkening coefficients, for which we assumed standard values (van Hamme 1993). Fourth, discretisation errors might be larger than expected, especially for the ETV, ECL data sets, where we used a linear interpolation between the neighboring time steps. We checked for all of the possibilities above, but none of them explain the tension.

3.3.3. Best fit.

How the best-fit solution with χ² = 3639180 looks is shown in Fig. 6. The fit of the SKY data set is acceptable; new WDS data of the fourth component are fitted perfectly, the NPOI data of the third component exhibit only minor systematics. Note that the coordinates are photometric with respect to Aa+Ab or Aa+Ab+B, respectively. The RVs are also acceptable; new CHIRON data are fitted perfectly, in agreement with line profiles (SYN); old RV data exhibit offsets up to 5 km/s, but not systematically. For the ETVs, the offsets are < 0.001 d, but one can notice a trend that MOST 2012 is fitted, TESS is fitted, but MOST 2017 is offset by −0.0015 d. It is clearly not easy to correct these issues, within the current dynamical model (see Sect. 3.4). We verified that the transversal light-time effect (Conroy et al. 2018; Vokrouhlický 2026) has a substantially smaller amplitude. Similar related offsets are seen in the ECL dataset, in agreement with light curves. The fits of VIS, CLO, and T3 data sets contain individual offsets for certain baselines, B/λ, but the trends of V²(B/λ) are correct. This is often caused by visibility calibration, or transfer-function issues. The synthetic SED data set is overestimated (by 0.1 to 0.15 mag), while SED2 is underestimated (by ∼0.3 mag). This might suggest an issue related to the fourth component, the total mass, m_sum, or the distance, d, but more likely to radiative parameters (see the discussion of T₃). The light-curve fit seems to be almost perfect for both the primary and secondary minima, but there are tiny, statistically significant offsets in minima timings, corresponding 1:1 to the ETV offsets above. The fit of normalised spectra (SYN) is also acceptable, including the Hα, Hβ wings, with a notable exception of the Hα, Hβ core depths for Aa, Ab components, which are slightly underestimated. This might suggest an issue related to the grids of synthetic spectra or relative luminosities of the components (Aa, Ab, B).

Fig. 6. Best-fit model from Fig. 5 with χ2 = 3639180 (unreduced). The light curves (LCs) were phased, to show different minima depths (MOST versus TESS). The synthetic spectra (SYN) plot presents only a subset of ten spectra and only the Hβ line region. One can still notice some systematics for ETV (MOST 2017 offset), ECL (TESS offset), VIS (individual observations), CLO (individual observations), SED (overestimated), SED2 (underestimated), or SYN (Aa, Ab line depths are not sufficient).

3.4. Models with tides and oblateness

In order to explain the trend seen in ETVs, in particular the negative offset −0.0015 d of the MOST 2017 timings (Fig. 6), we tested models with tides and oblateness². These terms induce a change of the period (Ṗ) and an additional precession ($\dot{\omega }$, $\dot{\mathrm{\Omega }}$), which is added to the existing precession due to four bodies. At the same time, we revised the uncertainties of ETVs to 10⁻⁴ d, in agreement with high-precision light curves and our aim to explain offsets at this level. We computed a grid of models parametrised by the oblateness, C₂₀, and the time lag, Δt (Mignard 1979). Only the closest components (Aa, Ab) are relevant, as the remaining (B, C) are too distant and the terms are negligible. For simplicity, we also assumed no misalignment and no spin precession.

The results in Fig. 7 show a clear signal for tides, but not for oblateness. The time lag, Δt, was about (800 ± 100) s when only the first orbit was free. If the second and third orbits were also free, the time lag could be as low as ∼100 s, but only at the expense of a poor fit of the ECL data set (Fig. 8).

Fig. 7. Same as Fig. 3 but for different fixed parameters, i.e., the oblateness, C20, 1, and the tidal time lag, Δt1. For the second component we assumed the same values. The uncertainties of ETVs were revised to 10−4 d. We tested 121 pairs of parameters, 300, 300, and 1000 sequential iterations of simplex: 1.9 × 105 models in total. The parameters of the inner Aa+Ab orbit (plus P2) were free. The best fit is indicated (red) along with a statistically equivalent, close-to-zero C20, 1 solution (green). Zero Δt1 is excluded.

Fig. 8. Same as Fig. 7 but the parameters of all orbits were free. The solutions differ from Fig. 7 partly because of perturbations by component C. The solutions with low Δt1 ∼ 100 s exhibit a poor fit of eclipse duration (ECL), so the solution with high Δt1 (and zero C20, 1) is still preferred.

We computed the tides for the reference radius R = R_⊙, the Love number k₂ = 0.3, the mean motion n ≐ 0.879 rad d⁻¹, and the rotation frequency ω = 1.25 rad d⁻¹. One can rescale Δt to any radius, for example R′=R₁ (and k₂′ = 0.03), because the tidal term is proportional to R⁵k₂Δt, so

$$\begin{array}{c}\hfill \mathrm{\Delta }t\prime =\mathrm{\Delta }t{(R/R\prime )}^5(k_2/k_2\prime )\simeq (710\pm 90)\mathrm{s}.\end{array}$$(6)

One can convert it to the dissipation factor

$$\begin{array}{c}\hfill Q=1/(\mathrm{\Delta }t2(\omega -n))\simeq 160,\end{array}$$(7)

or equivalently to Q′≡Q/k₂ ≃ 5500, or to $E\equiv \frac{2}{3}{k}_2\mathrm{\Delta }t2(\omega -n)\simeq 1.2\times {10}^{-4}$.

This seems to be too high for non-rotating B-type stars (Δt ≃ 0.15 s, Q ≃ 10⁶, Q′≃3 × 10⁷, E ≃ 2 × 10⁻⁸; Zahn 1975, 1977).

With this point of view in mind, the model with tides was excluded. We point out just a few arguments as to why dissipation could be larger than expected. Some stars orbited by exoplanets have substantially smaller Q′≃5 × 10⁴ (Penev et al. 2018; Maciejewski et al. 2016). Rotating stars exhibit Eddington–Sweet circulation (Eddington 1925; Sweet 1950), Soldberg–Høiland instability (Solberg 1936; Høiland 1941), or Goldreich–Schubert–Fricke instability (Goldreich & Schubert 1967; Fricke 1968), and such stars have smaller Q′. Some binaries show signs of large-amplitude, resonant oscillation, not only at the orbital frequency (‘heartbeat’), but at much higher multiples, for example 229 times (Fuller et al. 2017; Fuller & Felce 2024). Finally, the inner orbit (Aa+Ab) has non-zero, forced eccentricity, e → 0.008, due to component B, which increases dissipation as on Io (e = 0.004; Peale et al. 1979).

Alternative solutions also exist for non-zero oblateness, for example C_20, 1 = C_20, 2 = −0.001. Whether or not it is reasonable can be estimated from rotational and tidal (Roche) deformation. Assuming a Love number (Love 1909) of the order of k₂ ≃ 0.01 to 0.1 for radiative and convective stars, respectively (Claret 2004), the oblateness from rotation is (Kaula 1966)

$$\begin{array}{c}\hfill C_{20}=-\frac{2}{3}{k}_2{\left(\frac{\mathrm{\Omega }}{n}\right)}^2,\end{array}$$(8)

and from tides

$$\begin{array}{c}\hfill C_{20}=-\frac{2}{3}{k}_2{\left(\frac{R}{a}\right)}^3.\end{array}$$(9)

The maximal value (convective, rotational) is ∼1.5 × 10⁻⁴. This is why the model with oblateness was also excluded.

3.5. Models with five components ((Aa+Ab)+B)+(Ca+Cb))

Since some tension remains in the current model – specifically (i) high-precision astrometry (Tokovinin et al. 2020), leads to slightly more ‘elongated’ outer orbit; (ii) the RV dataset constrains the masses of components Aa, Ab, B; (iii) the ETV data set would require more massive component C; (iv) the SED dataset larger distance (69 pc); and (v) the SED2 data set a fainter component C – an elegant solution would be splitting of component C into a binary (Ca+Cb). This decreases its luminosity by up to ∼2.25 mag, depending on the mass ratio.

We computed the same grid of models as in Sect. 3.4. The results in Fig. 9 show a few solutions with non-zero tides, Δt = (100 ± 50) s, which is a more reasonable value than in Sect. 3.4, because it corresponds to Δt′ = 90 s, Q = 1300, Q′ = 44000, E = 1.5 × 10⁻⁵. Most importantly, though, this solution is no longer in contradiction with the ECL dataset, or the dynamics of the inner eclipsing binary. This confirms that a five-component model is necessary to correctly describe the ξ Tau system.

Fig. 9. Same as Fig. 7 but for models with five components (Aa+Ab+B+Ca+Cb). The perturbations and the light-time effect due to Ca+Cb components were stronger. The best-fit solution with low Δt1 ∼ 100 s now exhibits a good fit of eclipse duration (ECL).

4. Temporal evolution of ξ Tau

4.1. Signs of secular evolution

The ξ Tau system is appreciably compact to reveal temporal evolution of its individual orbits. In particular, the observations available to date require the orbital elements of the eclipsing pair (Aa+Ab) and of the third component (B) to evolve according to our dynamical model (Fig. 12). These variations occur over various timescales: (i) the shortest (days to a year) is defined by the orbital periods P₁ and P₂; (ii) the intermediate (decades to centuries) by the orbital period, P₃, or the precession rates ${\dot{\mathrm{\Omega }}}_1$, ${\dot{\mathrm{\Omega }}}_2$, ${\dot{\omega }}_2$, expressed in the Laplace reference frame, whose z-axis coincides with the angular momentum of the (Aa+Ab)+B triple; and (iii) the longest (millenia to 10⁵ years) is given by the slow precession rate ${\dot{\mathrm{\Omega }}}_3$.

4.1.1. Shortest timescales.

The mutual gravitational interaction in the (Aa+Ab)+B triple makes the shortest timescale apparent in the periods P₁ and P₂ (or equivalently, the semi-major axes a₁ and a₂) as well as the eccentricities e₁ and e₂. We find that e₁(t) oscillates between 0 and 0.008, which is dominated by short-period terms with P₁ and P₂ periods and amplitudes ≃0.003 and ≃0.004. These variations of e₁ are necessary to explain the exact eclipse timings, namely the contribution that cannot be produced by the perturbation of the binary’s mean motion by star B.

4.1.2. Intermediate timescales.

The correlated variations of inclinations i₁ and i₂, as well as the nodal longitudes Ω₁ and Ω₂ in the observer’s reference frame, occur with the periodicity of P_Ω ≃ 7000 d. The projected inclination, i₁, of the inner eclipsing binary changes from 86.1° to 87.1°, an effect clearly manifested in eclipse durations and depths, according to the precise photometric measurements by MOST (2012 and 2017) and TESS (from 2020 to 2023). These variations appear to be a consequence of a simple precession of the orbital angular momenta L_A and L_B of the Aa+Ab binary and of the component B orbit about their composite angular momentum, L_AB = L_A + L_B. In particular, both L_A and L_B describe coni with small opening angles of about L_AB with the period P_Ω, keeping their mutual angle J₁ nearly constant. Due to the very small value of e₁, the argument of periastron, ω₁, is an ill-defined orbital element (especially near epochs when e₁ ≃ 0), and its fast circulation is not a representative secular effect. In contrast, the argument of periastron, ω₂, of orbit B exhibits in the observer’s reference frame a regular precession with a rate of 2.1°  y⁻¹ (see Nemravová et al. 2016 and Fig. 12). This effect is evidenced by the measured RVs of components Aa, Ab, and B.

In order to explore the behaviour of the orbital eccentricity, e₁, of the eclipsing binary Aa+Ab on intermediate timescales, we digitally filtered the signal with periods P₁ and P₂ from the time series of osculating values of e₁. The secular theory of triple dynamics (e.g. Georgakarakos 2003, 2009; Breiter & Vokrouhlický 2015) predicts two contributions to e₁, namely (i) the free (proper) component, depending on the initial conditions, and (ii) the forced component, due to the presence of the third star (B) and its perturbations. We found that the forced component with an amplitude of 7.6 × 10⁻⁴ and a period of 167 yr (i.e. the apsidal precession timescale of ω₂) dominates the evolution of e₁, leaving no signal corresponding to the free component. This result suggests that the free component of e₁ has been damped by the past tidal evolution in the binary Aa+Ab, a process that is currently witnessed by the weak tidal signal (Sect. 3.4).

4.1.3. Longest timescales.

The variations over the longest timescales are associated with the outer orbit of component C (P₃ ≐ 18900 d ≐ 52 y). The gravitational perturbation of orbit B during the orbital motion of C is clearly manifested in the measured RVs. It is amplified near the periastron passage of C in July 2008 due to its large eccentricity, e₃ ≃ 0.573. An interesting ‘sawtooth’ structure seen during the sudden increase of i₂ between JD 2454000 and 2455000 (an interval of time extending approximately 1.5 y before and after the periastron passage of C; Fig. 12) may be explained as a coherent effect during five or six conjunctions between components B and C, repeating themselves after 145 d. We verified that the conjunctions occur at the same geometrical configuration, implying that the out-of-plane perturbing acceleration, W, in orbit B maintains the same sign. Consequently, since di/dt ∝ Wcos(f + ω) (Gauss 1809), the perturbing effect on i₂ accumulates coherently.

4.2. Stability

The gravitational perturbations also operate on timescales longer than centuries, beyond what is evidenced by the available observations. We thus conducted numerical integrations over 10⁵ y to study the stability of the ξ Tau system. We assumed the four-component model from Sect. 3.2. The integrator was the same, the Bulirsch–Stoer with an adaptive time step, but we used a smaller tolerance, ϵ = 10⁻¹⁰, to assure the angular momentum conservation.

We found that the ((Aa+Ab)+B)+C hierarchy is stable, with all components moving in a well-orchestrated manner. Apart from the (Aa+Ab)+B inner triple, where L_A and L_B precess about L_AB, as discussed in Sect. 4.1, we note that the nodal period P_Ω is smaller than P₃ (and non-resonant). This implies that the corresponding nodal period at the higher level in the ξ Tau hierarchy, namely the precession of L_AB and L_C about L_sum is much longer, i.e. ≃85000 y according to our model.

4.2.1. Non-existent Kozai oscillations.

For comparison, we tested a three-body model of ξ Tau, where the Aa+Ab binary was replaced with a single component A. This totally changed the overall orbital evolution, because of the Kozai oscillations (von Zeipel 1910; Lidov 1961; Kozai 1962). As illustrated in Fig. 11, the eccentricities e₂, e₃ and inclinations i₂, i₃ exhibit coupled oscillations, with a period of ≃11500 y. The peak value of e₂(t) is 0.8, which would destabilize the innermost binary Aa+Ab. Nevertheless, the stability of the (Aa+Ab)+B sub-system is ensured by a protection mechanism, consisting of its own mutual interaction, which induces the periastron precession, ${\dot{\omega }}_2$, with a period of only 180 y. This is much shorter than the Kozai timescale, effectively switching off the Kozai mechanism.

An interesting implication of the fact that L_C ≫ L_AB (see Fig. 10) is that the opening angle of the cone described by L_AB is relatively large (47°). As a result, the viewing geometry of the eclipsing binary (Aa+Ab) changes significantly over millennia. We predict eclipses will vanish in 18000 y from now.

Fig. 10. Angular-momentum vectors of the ξ Tau sub-systems (Aa+Ab, B, C) and their long-term evolution. The coordinates were rotated to the Laplace reference frame, with the z-axis along the total angular momentum, Lsum, of the whole ξ Tau system (in absence of external torques Lsum is constant). The individual contributions were obtained as ${\mathit{L}}_{\mathrm{A}}=m_{\mathrm{A}}\prime {\mathit{r}}_{\mathrm{A}}\times {\dot{\mathit{r}}}_{\mathrm{A}}$, where mA′ denotes the reduced mass, rA and ṙA the Jacobi coordinate and velocity of the Aa+Ab binary (and similarly for the B and C parts, LB and LC). The (2+1)+1 hierarchy of the system implies that (i) LC and LAB ≡ LA + LB regularly precess about Lsum with a period of ≃85000 y (see colours showing the time in the course of simulation), preserving their mutual angle J2 ≃ 71°; (ii) LA and LB regularly precess about LAB with a period of ≃7000 d, preserving their mutual angle of J1 ≃ 0.5°. Since LC holds the largest share of Lsum, its angular distance from Lsum is only 23°, while that of LAB is more than 47°. As a result, the observed inclinations with respect to the sky plane, i1 and i2, exhibit large variations over millennia.

Fig. 11. Same as Fig. 10 but only for three bodies (A, B, C), where A represents Aa+Ab. Without a protection mechanism, the Kozai oscillations were triggered, with a period of ≃11500 y. The angular momenta LB and LC still describe conical surfaces, but now their mutual angle J exhibits large oscillations, which triggers correlated oscillation of eccentricities and inclinations of the inner and outer orbits.

Fig. 12. Osculating orbital elements of ξ Tau in the course of time, for the model with χ2 = 3639180. The angular elements are expressed in the observer’s reference frame. See the description in the main text.

5. Discussion

5.1. Photocentre motion in Gaia DR4

In order to understand how the multiplicity of ξ Tau influences its parallax measurements (Vallenari et al. 2023), we computed its photocentre motion. The simplest, five-parameter model suitable for single stars consists of the position α, δ, the proper motion (PM) μ_α, μ_δ, and the parallax π. We added to this model the photocentre, computed as an average of barycentric positions, weighted by the passband fluxes of components Aa, Ab, B, and C. For reference, the Gaia G filter has λ_eff = 639.02 nm, Δ_eff = 317.32 nm and the respective weights were 0.181, 0.173, 0.603, and 0.042. Our results (Fig. 13) show that the overall photocentre motion within the observational time span of Gaia reaches almost 40 mas. However, this can be compensated by using different PM values, specifically μ_α′ = 67.0 mas y⁻¹, μ_δ′= − 40.5 mas y⁻¹. In other words, the new trajectory on the sky was adjusted so that its end points correspond to the old trajectory.

Fig. 13. Photocentre motion of the ξ Tau system and its influence on parallax measuments. We plot the orbital motion of individual components Aa, Ab, B, C in the barycentric frame (light gray or coloured, if within the observational timespan of Gaia), the photocentre motion (black), the parallactic motion (gray), the proper motion (light gray), and all contributions with the new parallax π′ = 14.717 mas (cyan). For comparison, we plot all contributions with the old parallax π = 16.791 mas (light cyan) and also the old trajectory from Gaia DR3 (dash-dotted).

Now the photocentre motion consists of two contributions: (i) 6 mas oscillating due to the triple (Aa+Ab)+B; and (ii) 4 mas curved due to component C. These offsets were most likely responsible for an incorrect value of π = 16.791 mas in Gaia DR3. When we used our preferred distance (67.9 pc), corresponding to π′ = 14.717 mas, the new trajectory became closer to the old trajectory, or the actual measurements. Consequently, our model explains not only all kinds of observations, but also why the Gaia DR3 parallax is offset. We do expect that binary, triple or quadruple astrometric models in Gaia DR4 will confirm our results.

5.2. Oscillations

In order to better understand rapid low-amplitude variability of ξ Tau (Nemravová et al. 2016), we used the new photometric data from TESS to compute a new periodogram. After excluding all eclipse intervals to suppress the orbital frequency, f_orb (and its harmonics), we actually computed a sequence of periodograms (Lenz & Breger 2005), revealing two dominant frequencies, f₁ = 1.407 d⁻¹, f₂ = 2.371 d⁻¹, with amplitudes of approximately 7 and 5 × 10⁻⁴ in the relative flux (Fig. 14). The corresponding periods are 0.710 and 0.421 d, respectively. The overall time span (1115 d) allowed us to reach the resolution of 1/Δ = 10⁻⁴ d⁻¹.

Fig. 14. Periodogram of ξ Tau from the TESS (top) and MOST 2012 (bottom) light curves. We removed eclipses and computed five Fourier transforms (Lenz & Breger 2005), with fitting of frequencies, amplitudes and phases, and subtracting the respective model, ∑iAisin[2π(fit + ϕi)]. For TESS, the dominant frequencies were 1.407, 2.371, 1.421, 2.814, and 0.279 d−1. We note f4 = 2f1 and f5 = 2forb of the inner orbit. For MOST 2012, on the contrary, we found 2.369, 2.216, 0.993, 0.282, and 2.839 d−1. The non-existent 1.407 d−1 frequency is labeled “ftess”.

The f₂ frequency is persistent and could be attributed to rotation of the third component (Nemravová et al. 2016). Its rotational broadening of v_rot3 = 242 km s⁻¹ corresponds to the upper limit of P_rot3 ≤ 0.472 d. The rotational axis inclination is unknown, but in order to obtain 0.421 d, it should be ∼60°.

The f₁ frequency should also be attributed to the third component, which is most likely a slowly pulsating B star (SPB; Waelkens & Rufener 1985). The star is located in the instability strip of g-modes of the order of ℓ = 1 or 2 (Sharma et al. 2022). However, there is a problem with non-existent f₁ oscillations in the MOST 2012 data, and no evidence of damping in the TESS 2020–2023 data. A possible explanation may consist of interference (‘beat’) of two nearby frequencies (e.g. Sharma et al. 2022, Fig. 5). However, the time resolution of the TESS photometry is insufficient to identify frequencies (f₁, f₁′) whose combination would result in conjectured beats (f₁ − f₁′) separated by more than eight years.

We verified that the reason for missing f₁ oscillations is not technical. The MOST telescope is a Maksutov–Cassegrain with a Fabry lens placed in the secondary focus (Walker et al. 2003). Its diaphragm (0.117 mm) corresponds to the angular diameter of ϕ = d/f = 27″, which is too large to miss light from the fourth component (up to 600 mas).

Still on a technical note, the TESS satellite observed ξ Tau as an over-exposed target, with blooming (Ricker et al. 2015; Krishnamurthy et al. 2019). In order to compute a sum of all pixels (both over- and under-exposed) the photometric aperture must be 13 pixels in the y-direction. Since the pixel size is 15 μm and the focal length is 146 mm, the angular distance is up to ϕ = 4.2′. Using Lightkurve (Cardoso et al. 2018), we checked stars in the surroundings. In particular, TIC 399947175 is close to the aperture, but it does not exhibit such oscillations. So, unfortunately, the f₁ oscillations remain an open problem.

6. Conclusions

In this work, we studied the ξ Tau quadruple system using dynamical models of differing complexities, constrained by 11 different types of observations. Using N-body and post-Newtonian terms and assuming four components (Aa, Ab, B, C), we found a global minimum of χ², where the total mass is 9.55 M_⊙ and the distance 66.1 pc (see Table 1). However, such a model still exhibits some tension; for example, minima timings for the MOST 2017 dataset are offset by ∼0.0015 d, the SED is overestimated by 30%, and so on.

It was surprisingly difficult to find a better solution. Eventually, we were forced to use N-body, post-Newtonian, oblateness, and tide terms and assume five components (Aa, Ab, B, Ca, Cb). We found a global minimum–where the total mass is 10.37 M_⊙–mainly due to the binary nature of component C and the distance 67.9 pc. Given the increased number of model parameters (36 vs. 44), it is not surprising that the model fits observations better. Nevertheless, the sole behavior of the model, how the minima timings were corrected, how the SED was corrected, suggests that the five-component model is the correct one.

Finally, let us point out that ξ Tau is indeed an interacting stellar system; even the most distant component C, interacts with component B at periastron passage, inducing a non-negligible phase shift, which propagates to the innermost eclipsing binary Aa+Ab, measurably shifting its eclipse timings. This has a potential to detect additional, dwarf, or exoplanetary components with low masses.

Acknowledgments

This work has been supported by the Czech Science Foundation through the grant 25-16507S (M. Brož and D. Vokrouhlický). We thank an anonymous referee for constructive comments. This work is co-funded by the European Union (ERC, MAGNIFY, Project 101126182). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work has made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC). This research has made use of the SIMBAD database operated at CDS, Strasbourg (France), and NASA’s Astrophysics Data System (ADS). This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA’s Science Mission directorate.

References

Appendix A: Rectification of CTIO/CHIRON spectra

As a prerequisite, one has to describe the CTIO/CHIRON instrument, in order to rectify all observed spectra and minimize systematics, in particular, in the Balmer lines wings. This required redefinition of the blaze function. The rectification is defined as

$$\begin{array}{c}\hfill I_i^{\mathrm{rect}}(\lambda )\equiv \frac{I_i^{\mathrm{obs}}(\lambda )}{R(\lambda )},\end{array}$$(A.1)

where I_i^rect(λ) is the rectified spectrum, I_i^obs(λ) the observed spectrum, R(λ) the rectification function; in the case of an echelle spectrum, it is called the ‘blaze function’. The motivation for rectification is to avoid absolute calibration, to measure line profiles, and to compare to normalized synthetic spectra. For complex echelle spectra, a number of methods exist: spline fitting, using a reference spectrum, a theoretical blaze function, α-shapes (Xu et al. 2019), or optimal extraction (Piskunov et al. 2021).

In the case of spline fitting (e.g. Harmanec et al. 2020), parts considered as continuum are fitted with a cubic spline, where the spline extends over wide lines. Among the common problems is uncertain continuum level, uncertain line wings, or unreliable log g values. Moreover, in the Paschen series region, overlapping lines and missing continuum make reliable rectification impossible.

In the case of a reference spectrum (Škoda et al. 2008), a rectified spectrum of a similar star is used. Alternatively, it could be a spectrum rectified by splines, obtained by comparison to a suitable synthetic spectrum. A measured spectrum is then divided by a reference spectrum (Eq. (A.1)), to get R(λ). Of course, this cannot be done without a low-noise reference spectrum. Among the common problems is noise in the reference spectrum, mismatch of resolution, or rapid changes.

Rectification by a theoretical blaze function (Barker 1984) assumes that

$$\begin{array}{c}\hfill R(\lambda )=A_r{\left\{\frac{sin\left[\pi {\alpha }_r{X}_r(\lambda )\right]}{\pi {\alpha }_r{X}_r(\lambda )}\right\}}^2\mathrm{for}\forall r,\end{array}$$(A.2)

where the coefficients are different for each order r, and

$$\begin{array}{c}\hfill X_r(\lambda )=r_{\mathrm{e}}(1-\frac{\lambda }{{\lambda }_{\mathrm{c},r}}),\end{array}$$(A.3)

where λ is the wavelength, λ_c the central wavelength, A the local maximum of the blaze function, r_e the echelle order; we use our own numbering of orders

$$\begin{array}{c}\hfill r\equiv 126-r_{\mathrm{e}}.\end{array}$$(A.4)

The grating constant

$$\begin{array}{c}\hfill \alpha =\frac{l}{d}cos\gamma ,\end{array}$$(A.5)

where l is the width of a grating facet, 1/d the spatial frequency of facets, γ the blaze angle.

Unfortunately, Eq. (A.2) only accounts for the basic grating theory. Slight variations are not included and it does not represent the actual CTIO/CHIRON data. The blaze function thus needs more degrees of freedom and our method is based on this approach.

A.1. Redefined blaze function

We proposed the following redefinition of the blaze function

$$\begin{array}{c}\hfill R(\lambda )=A_r{\left\{\frac{sin\left[\pi {\beta }_r(\lambda )X_r(\lambda )\right]}{\pi {\beta }_r(\lambda )X_r(\lambda )}\right\}}^2\mathrm{for}r=1..62,\end{array}$$(A.6)

where the constant α_r is replaced with the function

$$\begin{array}{cc}\hfill {\beta }_r(\lambda )& ={10}^{-10}{\gamma }_r{({\lambda }_{\mathrm{c},r}-{\epsilon }_r-\lambda )}^4+{10}^{-8}{\theta }_r{({\lambda }_{\mathrm{c},r}-{\epsilon }_r-\lambda )}^3+\hfill \\ \hfill & +{10}^{-6}{\delta }_r{({\lambda }_{\mathrm{c},r}-{\epsilon }_r-\lambda )}^2+{\alpha }_r,\hfill \end{array}$$(A.7)

while X_r(λ) remained the same as in Eq. (A.3); λ_c, r, A_r, α_r, δ_r, γ_r, θ_r, ε_r are order-dependent parameters, which we determined by fitting.

We defined the respective metric as

$$\begin{array}{c}\hfill {\chi }^2\equiv {\displaystyle \sum _{i=1}^N}{\left(\frac{I_i^{\mathrm{rect}}-I_i^{\mathrm{syn}}}{{\sigma }_i}\right)}^2,\end{array}$$(A.8)

where I_i^rect is the rectified observed intensity, I_i^syn is the intensity of the synthetic (rectified) spectrum; the uncertainties

$$\begin{array}{c}\hfill {\sigma }_i=\sqrt{\eta }\frac{I_i^{\mathrm{rect}}}{\sqrt{I_i^{\mathrm{obs}}}},\end{array}$$(A.9)

where I_i^obs is the observed (not rectified) intensity, η is the gain of the CCD detector.

Synthetic spectra were taken from our previous, converged model of ξ Tau (Nemravová et al. 2016). They served primarily to estimate the continuum level, not the depths of lines. We avoided the telluric lines, as obtained by the Skycalc tool (Noll et al. 2012; Jones et al. 2013) for the nearest available observatory, which is La Silla. The positions of the respective lines were masked in observed spectra.

In total, one has 2 parameters per-order and per-dataset (λ_c, r, A_r) and 5 parameters per-order (α_r, δ_r, γ_r, θ_r, ε_r), the number of orders is 62, the number of datasets is 10 (since we selected a subset, evenly distributed in time), which implied 1550 parameters. We optimized their values by the simplex algorithm (Nelder & Mead 1965), order by order. We call this first, direct method ‘resimplex1’. ³ The resulting blaze function is shown in Fig. A.1, and it is significantly better than our own, manual rectification.

Fig. A.1. Redefined blaze function R(λ), according to (Eq. (A.6)), estimated for CTIO/CHIRON spectra. It is plotted in absolute, analog-digital units (adu) and for shifted wavelengths λ − λc, r, for each order r. The rectification method was resimplex1. Colours correspond to χ2 contributions. The total χ2 = 666239 (unreduced) is a measure of correspondence between the rectified observed spectra and the synthetic (rectified) spectra.

A.2. Constrained blaze function

Alternatively, in order to decrease the variability of the blaze function, some parameters were approximated by polynomials

$$\begin{array}{cc}\hfill \alpha (r)& =a_1+a_{2}r+a_3{r}^2+a_4{r}^3+a_5{r}^4,\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill \gamma (r)& =g_1+g_{2}r+g_3{r}^2+g_4{r}^3+g_5{r}^4,\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill \delta (r)& =d_1+d_{2}r+d_3{r}^2+d_4{r}^3+d_5{r}^4,\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill \theta (r)& =t_1+t_{2}r+t_3{r}^2+t_4{r}^3+t_5{r}^4,\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill \epsilon (r)& =e_1+e_{2}r+e_3{r}^2+e_4{r}^3+e_5{r}^4,\hfill \end{array}$$(A.14)

while R(λ) and β_r(λ) functions remained the same. In the 1st step, the respective coefficients were estimated by fitting these polynomials to α_r, γ_r, δ_r, θ_r, ε_r discrete values (Fig. B.1). The remaining parameters λ_c, r, A_r were taken from Appendix A.1. In the 2nd step, we further optimized a₁, a₂, …, e₄, e₅, and λ_c, r, A_r were periodically updated during convergence, in order to improve their consistency with the polynomials. We again worked with 62 orders, 10 datasets, which implied 1265 parameters. We call this second, constrained method ‘resimplex2’.

In order to avoid problems with wide lines, we did not use the following orders for fitting (in Å): #7 4732.3–4786.8 Hβ, #8 4772.3–4827.3 Hβ, #9 4813.1–4868.5 Hβ, #10 4854.5–4910.5 Hβ, #11 4896.7–4953.2 Hβ, #39 6472.5–6547.0 Hα, #40 6547.7–6623.1 Hα, #44 6867.1–6946.1 telluric lines, #52 7609.5–7697.0 telluric lines, #57 8160.9–8254.7 telluric lines. The polynomials (Eqs. (A.10)–(A.14)) then serve as interpolants for these regions.

The result is shown in Fig. A.2. The χ² is higher, but still better than our manual rectification. The biggest contribution is from the last four orders (i.e. the Paschen series).

Fig. A.2. Same as Fig. A.1, but for the rectification method resimplex2. The respective χ2 = 891923.

A.3. Constrained central wavelengths

In order to further decrease the order-to-order variability, we constrained the central wavelength as

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{c}}(r)={\displaystyle \sum _{i=1}^{10}}{l}_i{r}^{i-1}\end{array}$$(A.15)

and also the amplitude as

$$\begin{array}{c}\hfill A(r)={\displaystyle \sum _{i=1}^9}{a}_i{r}^{i-1}.\end{array}$$(A.16)

This turned out to be the compromise between simplicity and complexity, keeping low χ². However, both λ_c(r) and A_r(r) must be different for each dataset, because they are directly related to the variability of the atmosphere (extinction, refraction). Other blaze parameters were taken from Appendix A.2. We call this third method ‘resimplex3’. For technical reasons, spectra had to be split at orders 16 (CCD amplifiers), 28, 48 (Paschen series). We processed all 277 datasets this way, but we worked with the datasets sequentially, which limited the number of parameters.

The resulting blaze function is shown in Fig. A.3. Specifically, how the rectification works in the Hα line region is shown in Fig. 1. The respective χ² value is even higher, but still acceptable, if one avoids the Paschen series region. In fact, the relatively higher χ² value does not imply worse rectification. Instead, interpolating over some orders (Hα, Hβ) implies lower systematic errors due to the uncertain continuum level.

Fig. A.3. Same as Fig. A.1, but for the rectification method resimplex3. The respective χ2 = 1217823.

Appendix B: Supplementary figures

Here we show parameters of the blaze function (Fig. B.1), synthetic spectra of individual components Aa, Ab, B, C (Fig. B.2), the ETVs of ξ Tau for the model with tides (Fig. B.3), and the MCMC analysis for the model without tides (Fig. B.4).

Fig. B.1. Parameters of the blaze function as a function of the order r, shown either as discrete values αr, δr, ϵr, θr, γr, or as continuous polynomials α(r), δ(r), ϵ(r), θ(r), γ(r). We show a comparison of two rectification methods: resimplex1 on a set of 10 spectra (black), resimplex1 on a different set of 10 spectra (gray), resimplex2 (blue), resimplex2 with splitting at orders 16, 28, 48 (red). We also show a preliminary fit of black values (green), which served as a starting for resimplex2. The variability of discrete values is substantial, especially for orders r ≳ 50 (Paschen series), but this is suppressed when using polynomials.

Fig. B.2. Synthetic spectra of individual components Aa, Ab, B, C for the model with χ2 = 3639180, in the Hβ line region.

Fig. B.3. Same as Fig. 2 for non-zero tides, with revised uncertainties of ETVs of the order of 10−4 d. The total χ2 = 4145. Systematics of ETVs were largely compensated by suitable a value of Δt.

Fig. B.4. MCMC analysis for the reference model with χ2 = 1553, constrained by SKY, RV, ETV, ECL, SED, SED2 datasets. The corner plot shows uncertainties and correlations of 30 parameters, which were sampled with 64 walkers, 1000 iterations, with 500 burn-in steps. The highest-probability model is plotted (green). Specifically, msum = 9.572 ± 0.023 M⊙, q1 = 0.949 ± 0.001, q2 = 0.851 ± 0.004, q3 = 0.176 ± 0.001, P1 = 7.14730 ± 0.00010 d, log e1 = −2.835 ± 0.014, i1 = 87.050 ± 0.204°, Ω1 = 328.727 ± 0.143°, ϖ1 = 113.833 ± 0.646°, λ1 = 58.676 ± 0.144°, P2 = 145.779 ± 0.005 d, log e2 = −0.682 ± 0.003, i2 = 86.753 ± 0.120°, Ω2 = 328.575 ± 0.133°, ϖ2 = 340.150 ± 0.666°, λ2 = 63.402 ± 0.223°, P3 = 18803.7 ± 68.1 d, log e3 = −0.242 ± 0.001, i3 = −23.899 ± 0.157°, Ω3 = 101.493 ± 0.642°, ϖ3 = 117.655 ± 0.196°, λ3 = 148.337 ± 0.088°, T1 = 10520 ± 39 K, T2 = 10369 ± 77 K, T3 = 14060 ± 107 K, log g1 = 4.366 ± 0.019, log g2 = 4.352 ± 0.015, log g3 = 4.302 ± 0.008, γ = 10.132 ± 0.070 km s−1, d = 65.951 ± 0.112 pc. These 1-σ uncertainties are considered local and small, compared to global ones in Table 1, because of some systematics and tension remained in the model.

Fig. B.5. Comparison of CTIO/CHIRON observed spectra (one per night) with synthetic spectra, for the model with χ2 = 3639180. The fit is acceptable, including the Hα, Hβ wings, with a notable exception of the Hα, Hβ core depths for Aa, Ab components (see red circles for spectra with higher signal-to-noise ratio).

Appendix C: Supplementary tables

Here we show the journal of photometric observations (Table C.1) and dependent parameters of ξ Tau models (Table C.2)

All Tables

All Figures

Fig. 1. Rectification of a CTIO/CHIRON echelle spectrum (ktc00027) in the region of the Hα line (6563 Å), using the resimplex3 method. Top: Synthetic spectra were computed from our model in Table 1, ‘all-data’ column. Four component spectra (Aa, Ab, B, C) and the resulting spectrum are plotted. Second row: Observed spectrum from CHIRON (black), orders r = 39 to 41, and the corresponding blaze function (green). Third row: Rectified spectrum (black), its comparison to the synthetic spectrum (blue), and corresponding residuals (red). The rectification was done with respect to the continuum level, with interpolation of the blaze parameters over wide lines (cyan) and problematic orders, also avoiding locations of the telluric lines (green). Bottom: Residuals plotted separately.

In the text

Fig. 2. Reference model of ξ Tau, showing ETVs (top), eclipse durations (middle), perturbations with respect to a two-body ephemeris (bottom), MOST and TESS observations (blue), model (grey), residuals (red), and, for reference, observed minima (green). The total χ2 = 1553 (see Table 1 for individual contributions). The mean ephemeris used in the bottom panel is P = 7.146665 d, T0 = 2456224.704705. Our model is certainly able to explain the major perturbations, which are substantial (±0.02 d, or ±30 min), but some systematics remain to be explained (up to 0.002 d, or 3 min).

In the text

Fig. 3. Best-fit models of ξ Tau constrained by multi-technique observations for two fixed parameters, the total mass msum and the distance d. The χ2 values (colours) indicate best fit (cyan), good fits (blue), and poor fits (orange); the overall best fit (red) has χ2 = 1661 (unreduced). We tested 121 pairs of parameters and 3000 sequential iterations of simplex, i.e. 3.6 × 105 models in total. Other dynamical parameters (except msum, d) were free; other radiative parameters were fixed (Nemravová et al. 2016). We constrained the model by SKY, RV, ETV, ECL, SED, SED2 data sets; we did not use interferometry (but we used SKY), spectroscopy (but RV), or light curves (but ETV, ECL). We assumed unit weights, except wetv = wecl = 100, because we wanted to exclude all models with incorrect dynamics. The two parameters are often positively correlated (e.g. SKY, VIS, CLO, T3), as predicted (Kepler 1619), but some data sets are distance independent or orthogonal (RV, ETV, SED, SED2), and the model is considered well constrained. After re-convergence (green), the χ2 further decreased to 1553. The Gaia distance of 59.6 pc was excluded. Additionally, VIS, CLO, and T3 were also computed, but their weights were set to zero, otherwise our model would exhibit tension between VIS and CLO.

In the text

	Fig. 4. Best-fit model from Fig. 3 with χ2 = 1553 (unreduced). Plots for each data set, SKY, RV, ETV, ECL, VIS, CLO, T3, SED, and SED2, contain observations (blue), synthetic data (yellow), and residuals (red). The best-fit distance is d = 66.0 pc. One can notice some systematics for ETV (see MOST 2017), VIS, or CLO data sets.
In the text

	Fig. 5. Same as Fig. 3 but for two different fixed parameters, the effective temperatures T1and T3. This time, light curves (LC) and synthetic spectra (SYN) were included. The overall best fit χ2 = 3639180 (unreduced). One can notice some systematics; for example, there is tension between VIS and CLO, the best ECL is offset, the best SED is offset, and the best SED2 is very offset.
In the text

Fig. 6. Best-fit model from Fig. 5 with χ2 = 3639180 (unreduced). The light curves (LCs) were phased, to show different minima depths (MOST versus TESS). The synthetic spectra (SYN) plot presents only a subset of ten spectra and only the Hβ line region. One can still notice some systematics for ETV (MOST 2017 offset), ECL (TESS offset), VIS (individual observations), CLO (individual observations), SED (overestimated), SED2 (underestimated), or SYN (Aa, Ab line depths are not sufficient).

In the text

Fig. 7. Same as Fig. 3 but for different fixed parameters, i.e., the oblateness, C20, 1, and the tidal time lag, Δt1. For the second component we assumed the same values. The uncertainties of ETVs were revised to 10−4 d. We tested 121 pairs of parameters, 300, 300, and 1000 sequential iterations of simplex: 1.9 × 105 models in total. The parameters of the inner Aa+Ab orbit (plus P2) were free. The best fit is indicated (red) along with a statistically equivalent, close-to-zero C20, 1 solution (green). Zero Δt1 is excluded.

In the text

	Fig. 8. Same as Fig. 7 but the parameters of all orbits were free. The solutions differ from Fig. 7 partly because of perturbations by component C. The solutions with low Δt1 ∼ 100 s exhibit a poor fit of eclipse duration (ECL), so the solution with high Δt1 (and zero C20, 1) is still preferred.
In the text

	Fig. 9. Same as Fig. 7 but for models with five components (Aa+Ab+B+Ca+Cb). The perturbations and the light-time effect due to Ca+Cb components were stronger. The best-fit solution with low Δt1 ∼ 100 s now exhibits a good fit of eclipse duration (ECL).
In the text

Fig. 10. Angular-momentum vectors of the ξ Tau sub-systems (Aa+Ab, B, C) and their long-term evolution. The coordinates were rotated to the Laplace reference frame, with the z-axis along the total angular momentum, Lsum, of the whole ξ Tau system (in absence of external torques Lsum is constant). The individual contributions were obtained as ${\mathit{L}}_{\mathrm{A}}=m_{\mathrm{A}}\prime {\mathit{r}}_{\mathrm{A}}\times {\dot{\mathit{r}}}_{\mathrm{A}}$, where mA′ denotes the reduced mass, rA and ṙA the Jacobi coordinate and velocity of the Aa+Ab binary (and similarly for the B and C parts, LB and LC). The (2+1)+1 hierarchy of the system implies that (i) LC and LAB ≡ LA + LB regularly precess about Lsum with a period of ≃85000 y (see colours showing the time in the course of simulation), preserving their mutual angle J2 ≃ 71°; (ii) LA and LB regularly precess about LAB with a period of ≃7000 d, preserving their mutual angle of J1 ≃ 0.5°. Since LC holds the largest share of Lsum, its angular distance from Lsum is only 23°, while that of LAB is more than 47°. As a result, the observed inclinations with respect to the sky plane, i1 and i2, exhibit large variations over millennia.

In the text

Fig. 11. Same as Fig. 10 but only for three bodies (A, B, C), where A represents Aa+Ab. Without a protection mechanism, the Kozai oscillations were triggered, with a period of ≃11500 y. The angular momenta LB and LC still describe conical surfaces, but now their mutual angle J exhibits large oscillations, which triggers correlated oscillation of eccentricities and inclinations of the inner and outer orbits.

In the text

	Fig. 12. Osculating orbital elements of ξ Tau in the course of time, for the model with χ2 = 3639180. The angular elements are expressed in the observer’s reference frame. See the description in the main text.
In the text

Fig. 13. Photocentre motion of the ξ Tau system and its influence on parallax measuments. We plot the orbital motion of individual components Aa, Ab, B, C in the barycentric frame (light gray or coloured, if within the observational timespan of Gaia), the photocentre motion (black), the parallactic motion (gray), the proper motion (light gray), and all contributions with the new parallax π′ = 14.717 mas (cyan). For comparison, we plot all contributions with the old parallax π = 16.791 mas (light cyan) and also the old trajectory from Gaia DR3 (dash-dotted).

In the text

Fig. 14. Periodogram of ξ Tau from the TESS (top) and MOST 2012 (bottom) light curves. We removed eclipses and computed five Fourier transforms (Lenz & Breger 2005), with fitting of frequencies, amplitudes and phases, and subtracting the respective model, ∑iAisin[2π(fit + ϕi)]. For TESS, the dominant frequencies were 1.407, 2.371, 1.421, 2.814, and 0.279 d−1. We note f4 = 2f1 and f5 = 2forb of the inner orbit. For MOST 2012, on the contrary, we found 2.369, 2.216, 0.993, 0.282, and 2.839 d−1. The non-existent 1.407 d−1 frequency is labeled “ftess”.

In the text

Fig. A.1. Redefined blaze function R(λ), according to (Eq. (A.6)), estimated for CTIO/CHIRON spectra. It is plotted in absolute, analog-digital units (adu) and for shifted wavelengths λ − λc, r, for each order r. The rectification method was resimplex1. Colours correspond to χ2 contributions. The total χ2 = 666239 (unreduced) is a measure of correspondence between the rectified observed spectra and the synthetic (rectified) spectra.

In the text

	Fig. A.2. Same as Fig. A.1, but for the rectification method resimplex2. The respective χ2 = 891923.
In the text

	Fig. A.3. Same as Fig. A.1, but for the rectification method resimplex3. The respective χ2 = 1217823.
In the text

Fig. B.1. Parameters of the blaze function as a function of the order r, shown either as discrete values αr, δr, ϵr, θr, γr, or as continuous polynomials α(r), δ(r), ϵ(r), θ(r), γ(r). We show a comparison of two rectification methods: resimplex1 on a set of 10 spectra (black), resimplex1 on a different set of 10 spectra (gray), resimplex2 (blue), resimplex2 with splitting at orders 16, 28, 48 (red). We also show a preliminary fit of black values (green), which served as a starting for resimplex2. The variability of discrete values is substantial, especially for orders r ≳ 50 (Paschen series), but this is suppressed when using polynomials.

In the text

	Fig. B.2. Synthetic spectra of individual components Aa, Ab, B, C for the model with χ2 = 3639180, in the Hβ line region.
In the text

	Fig. B.3. Same as Fig. 2 for non-zero tides, with revised uncertainties of ETVs of the order of 10−4 d. The total χ2 = 4145. Systematics of ETVs were largely compensated by suitable a value of Δt.
In the text

Fig. B.4. MCMC analysis for the reference model with χ2 = 1553, constrained by SKY, RV, ETV, ECL, SED, SED2 datasets. The corner plot shows uncertainties and correlations of 30 parameters, which were sampled with 64 walkers, 1000 iterations, with 500 burn-in steps. The highest-probability model is plotted (green). Specifically, msum = 9.572 ± 0.023 M⊙, q1 = 0.949 ± 0.001, q2 = 0.851 ± 0.004, q3 = 0.176 ± 0.001, P1 = 7.14730 ± 0.00010 d, log e1 = −2.835 ± 0.014, i1 = 87.050 ± 0.204°, Ω1 = 328.727 ± 0.143°, ϖ1 = 113.833 ± 0.646°, λ1 = 58.676 ± 0.144°, P2 = 145.779 ± 0.005 d, log e2 = −0.682 ± 0.003, i2 = 86.753 ± 0.120°, Ω2 = 328.575 ± 0.133°, ϖ2 = 340.150 ± 0.666°, λ2 = 63.402 ± 0.223°, P3 = 18803.7 ± 68.1 d, log e3 = −0.242 ± 0.001, i3 = −23.899 ± 0.157°, Ω3 = 101.493 ± 0.642°, ϖ3 = 117.655 ± 0.196°, λ3 = 148.337 ± 0.088°, T1 = 10520 ± 39 K, T2 = 10369 ± 77 K, T3 = 14060 ± 107 K, log g1 = 4.366 ± 0.019, log g2 = 4.352 ± 0.015, log g3 = 4.302 ± 0.008, γ = 10.132 ± 0.070 km s−1, d = 65.951 ± 0.112 pc. These 1-σ uncertainties are considered local and small, compared to global ones in Table 1, because of some systematics and tension remained in the model.

In the text

	Fig. B.5. Comparison of CTIO/CHIRON observed spectra (one per night) with synthetic spectra, for the model with χ2 = 3639180. The fit is acceptable, including the Hα, Hβ wings, with a notable exception of the Hα, Hβ core depths for Aa, Ab components (see red circles for spectra with higher signal-to-noise ratio).
In the text