An advanced N-body model for interacting multiple stellar systems

TL;DR

This paper introduces an advanced N-body simulation model for complex stellar systems that integrates multiple observational data types to produce comprehensive and artifact-free models, enabling detailed dynamical studies.

Contribution

The paper presents a novel integrated N-body modeling approach combining various observational data and a modified Wilson--Devinney code for comprehensive analysis of multiple stellar systems.

Findings

Robust models can be constructed using joint observational data and optimization.

The N-body model avoids artifacts from non-self-consistent gravitational interactions.

Examples demonstrate the model's capability to study dynamical effects.

Abstract

We construct an advanced model for interacting multiple stellar systems in which we compute all trajectories with a numerical N-body integrator, namely the Bulirsch--Stoer from the SWIFT package. We can then derive various observables: astrometric positions, radial velocities, minima timings (TTVs), eclipse durations, interferometric visibilities, closure phases, synthetic spectra, spectral-energy distribution, and even complete light curves. We use a modified version of the Wilson--Devinney code for the latter, in which the instantaneous true phase and inclination of the eclipsing binary are governed by the N-body integration. If one has all kinds of observations at disposal, a joint $\chi^2$ metric and an optimisation algorithm (a~simplex or simulated annealing) allows to search for a global minimum and construct very robust models of stellar systems. At the same time, our N-body model is free from artefacts which may arise if mutual gravitational interactions among all components are not self-consistently accounted for. Finally, we present a number of examples showing dynamical effects that can be studied with our code and we discuss how systematic errors may affect the results (and how to prevent this from happening).

Tables

Table 1. Table 1: Notation used for coordinates, velocities, and a number of other quantities and uncertainties, which we use in our N-body model.

$N_{\mathrm{bod}}$	number of bodies
$m$	mass ($G{M}_{\odot }$ units)
$q=\frac{m_1}{m_2}$	mass ratio
${\eta }_{ij}=\frac{m_j{m}_i}{{\left(m_j+m_i\right)}^2}$	symmetrized mass ratio
$k_L$	Love number
${\omega }_{\mathrm{rot}}$	rotational angular velocity
$x_b,y_b,z_b$	barycentric coordinates
$v_{xb},v_{yb},v_{zb}$	barycentric velocities
$x_h,y_h,z_h$	1-centric coordinates
$v_{xh},v_{yh},v_{zh}$	1-centric velocities
$x_p,y_p$	1+2 photocentric sky-plane coordinates
$x_{\mathrm{p3}},y_{\mathrm{p3}}$	1+2+3 photocentric coordinates
$x_a=\frac{x_h}{d},y_a$	1-centric coordinates in an angular measure
$\hat{X},\hat{Y},\hat{Z}$	unitvectors aligned with 1+2 eclipsing pair
$\hat{O}=(0,0,-1)$	observers direction
$\gamma$	systemic velocity
$v_{\mathrm{rad}}$	observed radial velocity
$t_{\mathrm{ecl}}$	mid-epoch of an eclipse of 1+2 pair
${ϵ}_{\mathrm{ecl}}$	eclipse duration
$L,L_{\mathrm{tot}}$	component luminosity and the total one
$T_{\mathrm{eff}}$	effective temperature
$R$	stellar radius

Table 2. Table 1: Cont.

$\lambda$, $\Delta \lambda$	effective wavelength and bandwidth
$B_{\lambda }\left(T\right)$	the Planck function
$V$	complex visibility, squared visibility is ${\left|V\right|}^2$
$T_3$	complex triple product, closure phase is $\arg T_3$
$u,v$	projected baselines (expressed in cycles, $\frac{B}{\lambda }$)
$\theta =\frac{2R}{d}$	angular diameter
$u_{\mathrm{limb}}$	linear limb-darkening coefficient
$d$	distance to the system
$m_V$	magnitude (in V band or another)
$m_0$	zero point
$I_{\lambda }$, $I_{\mathrm{syn}}$	normalized monochromatic intensity
$F_{\mathrm{syn}}$	absolute monochromatic flux (in $\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-2}{\mathrm{cm}}^{-1}$)
$F_{V\mathrm{calib}}$	calibration flux
$f_V$	filter transmission coefficient
$g=\frac{GM}{R^2}$	surface gravity, $\log g$ in cgs
$v_{\mathrm{rot}}$	projected rotational velocity
$𝒵$	metallicity
${\sigma }_{\mathrm{sky}\mathrm{major},\mathrm{minor}}$	uncertainty of the astrometric position,
	angular sizes of the uncertainty ellipse
${\varphi }_{\mathrm{ellipse}}$	position angle of the ellipse
$𝐑\left(\dots \right)$	the corresponding $2\times 2$ rotation matrix
${\sigma }_{\mathrm{rv}}$	uncertainty of the radial velocity
${\sigma }_{\mathrm{ttv}}$	uncertainty of the eclipse mid-epoch timing

Table 3. Table 1: Cont.

${\sigma }_{\mathrm{ecl}}$	uncertainty of the eclipse duration
${\sigma }_{\mathrm{vis}}$	uncertainty of the squared visibility
${\sigma }_{\mathrm{clo}}$	uncertainty of the closure phase
${\sigma }_{\mathrm{t3}}$	uncertainty of the triple product amplitude
${\sigma }_{\mathrm{lc}}$	uncertainty of the light-curve data
${\sigma }_{\mathrm{syn}}$	uncertainty of the normalized intensity
${\sigma }_{\mathrm{sed}}$	uncertainty of the spectral-energy distribution
$m_j^{\min },m_j^{\max }$	minimum and maximum masses

Table 4. Table 2: The mock system parameters which were used to generate synthetic data for testing. The notation is the same as in Table  1 . Values rounded to typical uncertainties of the parameters are presented in this table. The osculating elements correspond to the epoch T 0 = 2456224.724705 subscript 𝑇 0 2456224.724705 T_{0}=2456224.724705 .

Par.	Value							Unit
$m_1$	$2.238483$	$m_2$	$2.009645$	$m_3$	$3.7472$	$m_4$	$0.92$	$M_{\odot }$
$a_1$	$0.1176356$	$a_2$	$1.08420$	$a_3$	$28.39$			au
$e_1$	$0.0000$	$e_2$	$0.2167$	$e_3$	$0.568$
$i_1$	$87.6$	$i_2$	$86.3$	$i_3$	$-18.2$			deg
${\Omega }_1$	$329.1$	${\Omega }_2$	$328.6$	${\Omega }_3$	$114.7$			deg
${\omega }_1$	$275.65$	${\omega }_2$	$0.00$	${\omega }_3$	$1.0$			deg
$M_1$	$174.44$	$M_2$	$88.04$	$M_3$	$32.7$			deg
$\gamma$	$8.82$							$\mathrm{km}{s}^{-1}$
$d$	$67.9$							pc
$T_{\mathrm{eff1}}$	$10727$	$T_{\mathrm{eff2}}$	$10275$	$T_{\mathrm{eff3}}$	$13120$	$T_{\mathrm{eff4}}$	$6526$	K
$R_1$	$1.586$	$R_2$	$1.642$	$R_3$	$2.727$	$R_4$	$0.877$	$R_{\odot }$
$v_{\mathrm{rot1}}$	$16.2$	$v_{\mathrm{rot2}}$	$12.5$	$v_{\mathrm{rot3}}$	$234.9$	$v_{\mathrm{rot4}}$	$80.1$	$\mathrm{km}{s}^{-1}$
$m_{01}$	$1.000$	$m_{02}$	$3.335$	$m_{03}$	$3.646$	$m_{04}$	$3.730$	mag

Table 5. Table 3: Suggested examples of compact stellar systems for which the N-body model could be useful (or inevitable). This ‘catalogue’ obviously cannot be considered comprehensive.

Designation	Reference
$\lambda$ Tau	Fekel & Tomkin (1982)
$\xi$ Tau	Nemravová et al. (2016)
VW LMi	Pribulla et al. (2008)
V994 Her = HD 170314	Zasche & Uhlář (2016)
V907 Sco = HD 163302	Lacy et al. (1999)
HD 91962	Tokovinin et al. (2015)
HD 109648	Jha et al. (2000)
HD 144548	Alonso et al. (2015)
HD 181068 = KIC 5952403	Fuller et al. (2013)
KIC 05255552	Borkovits et al. (2016)
KIC 05771589
KIC 06964043
KIC 07289157
KIC 07668648
KIC 07955301
KIC 09714358

Equations

$${\ddot{\mathbf{r}}}_{\mathrm{b}i}=-\sum _{j\ne i}^{N_{\mathrm{bod}}}\frac{G{m}_j}{r_{ji}^3}{\mathbf{r}}_{ji}+{\mathbf{f}}_{\mathrm{tidal}}+{\mathbf{f}}_{\mathrm{oblat}}+{\mathbf{f}}_{\mathrm{ppn}},$$

$${\mathbf{f}}_{\mathrm{tidal}}=-\sum _{j\ne i}^{N_{\mathrm{bod}}}3k_{\mathrm{L}i}\frac{G{m}_j^2}{m_i}\frac{R_i^5}{r_{ij}^8}{\mathbf{r}}_{ji},$$

$${\mathbf{f}}_{\mathrm{oblat}}=-\sum _{j\ne i}^{N_{\mathrm{bod}}}\frac{1}{2}{k}_{\mathrm{L}j}{\omega }_{\mathrm{rot}j}^2\frac{R_i^5}{r_{ij}^5}{\mathbf{r}}_{ji},$$

$${\mathbf{f}}_{\mathrm{ppn}}$$

$$t_{\mathrm{ecl}}^{\prime }=t_{\min }+\frac{z_{\mathrm{b}1+2}-z_{\mathrm{b}1+2}(t=T_0)}{c}-\frac{z_{\mathrm{h}2}}{c},$$

$${ϵ}_{\mathrm{ecl}}^{\prime }=\frac{2}{{\bar{v}}_{\mathrm{h}2}}\sqrt{(R_1+R_2{)}^2-{\Delta }_{\min }^2},$$

$$L_j(T_{\mathrm{eff}j},R_j)=4\pi R_j^2{\int }_{\lambda -\Delta \lambda /2}^{\lambda +\Delta \lambda /2}{F}_{\mathrm{syn}}\left(\lambda ,T_{\mathrm{eff}j},\log g_j,v_{\mathrm{rot}j},{\mathcal{Z}}_j\right)\mathrm{d}\lambda ,$$

$$V^{\prime }(u,v)$$

$$T_3^{\prime }=V^{\prime }(u_1,v_1)V^{\prime }(u_2,v_2)V^{\prime }(-(u_1+u_2),-(v_1+v_2)),$$

$${\phi }_{\mathrm{ecl}}^{\prime }=\frac{1}{2\pi }\arctan \frac{\hat{O}\cdot \hat{Y}}{\hat{O}\cdot \hat{X}},$$

$$i_{\mathrm{ecl}}^{\prime }=\arccos (-\hat{O}\cdot \hat{Z}),$$

$${\Omega }_{\mathrm{Kopal}j}\doteq {⟨\frac{1}{r_1}+\frac{q}{r_2}+\frac{1}{2}(1+q)r_3^2⟩}_{\mathrm{circle}\frac{R_j}{r_{12}},\mathrm{where}}$$

$${\mathbf{r}}_1\in \mathrm{circle},{\mathbf{r}}_2={\mathbf{r}}_1-(1,0,0),{\mathbf{r}}_3=(x_1-\frac{q}{1+q},y_1,0),$$

$$I_{\lambda }^{\prime }=\sum _{j=1}^{N_{\mathrm{bod}}}\frac{L_j}{L_{\mathrm{tot}}}{I}_{\mathrm{syn}}\left[\lambda \left(1-\frac{v_{z\mathrm{b}j+\gamma }}{c}\right),T_{\mathrm{eff}j},\log g_j,v_{\mathrm{rot}j},{\mathcal{Z}}_j\right],$$

$$F_V^{\prime }=\sum _{j=1}^{N_{\mathrm{bod}}}{\left(\frac{R_j}{d}\right)}^2{\int }_0^{\infty }{F}_{\mathrm{syn}}\left[\lambda ,T_{\mathrm{eff}j},\log g_j,v_{\mathrm{rot}j},{\mathcal{Z}}_j\right]f_V(\lambda )\mathrm{d}\lambda ,$$

$$m_V^{\prime }=-2.5{\log }_{10}\frac{F_V^{\prime }}{F_{V\mathrm{calib}}{\int }_0^{\infty }{f}_V(\lambda )\mathrm{d}\lambda },$$

$${\chi }^2$$

$${\chi }_{\mathrm{sky}}^2=\sum _{j=1}^{N_{\mathrm{bod}}}\sum _{i=1}^{N_{\mathrm{sky}j}}\left\{\frac{(\Delta x_{ji}{)}^2}{{\sigma }_{\mathrm{sky}\mathrm{major}ji}^2}+\frac{(\Delta y_{ji}{)}^2}{{\sigma }_{\mathrm{sky}\mathrm{minor}ji}^2}\right\},$$

$$(\Delta x_{ji},\Delta y_{ji})=\mathbf{R}\left(-{\varphi }_{\mathrm{ellipse}}-\frac{\pi }{2}\right)\times \text{pmatrix}{x}_{\mathrm{p}ji}^{\prime }-x_{\mathrm{p}ji}\text{cr}{y}_{\mathrm{p}ji}^{\prime }-y_{\mathrm{p}ji},$$

$${\chi }_{\mathrm{rv}}^2=\sum _{j=1}^{N_{\mathrm{bod}}}\sum _{i=1}^{N_{\mathrm{rv}j}}\frac{{\left(v_{z\mathrm{b}ji}^{\prime }+\gamma -v_{\mathrm{rad}ji}\right)}^2}{{\sigma }_{\mathrm{rv}ji}^2},$$

$${\chi }_{\mathrm{ttv}}^2=\sum _{i=1}^{N_{\mathrm{ttv}}}\frac{{\left(t_{\mathrm{ecl}i}^{\prime }-t_{\mathrm{ecl}i}\right)}^2}{{\sigma }_{\mathrm{ttv}i}^2},$$

$${\chi }_{\mathrm{ecl}}^2=\sum _{i=1}^{N_{\mathrm{ecl}}}\frac{{\left({ϵ}_{\mathrm{ecl}i}^{\prime }-{ϵ}_{\mathrm{ecl}i}\right)}^2}{{\sigma }_{\mathrm{ecl}i}^2},$$

$${\chi }_{\mathrm{vis}}^2=\sum _{i=1}^{N_{\mathrm{vis}}}\frac{{\left({|V^{\prime }(u_i,v_i)|}^2-|V{|}_i^2\right)}^2}{{\sigma }_{\mathrm{vis}i}^2},$$

$${\chi }_{\mathrm{clo}}^2=\sum _{i=1}^{N_{\mathrm{clo}}}\frac{{\left(\arg T_{3i}^{\prime }-\arg T_{3i}\right)}^2}{{\sigma }_{\mathrm{clo}i}^2},$$

$${\chi }_{t3}^2=\sum _{i=1}^{N_{t3}}\frac{{\left(|T_{3i}{|}^{\prime }-|T_{3i}|\right)}^2}{{\sigma }_{t3i}^2},$$

$${\chi }_{\mathrm{lc}}^2=\sum _{k=1}^{N_{\mathrm{band}}}\sum _{i=1}^{N_{\mathrm{lc}k}}\frac{{\left(m_{Vki}^{\prime }+m_{0k}-m_{Vki}\right)}^2}{{\sigma }_{\mathrm{lc}ki}^2},$$

$${\chi }_{\mathrm{syn}}^2=\sum _{i=1}^{N_{\mathrm{syn}}}\frac{{\left(I_{\lambda i}^{\prime }-I_{\lambda i}\right)}^2}{{\sigma }_{\mathrm{syn}i}^2},$$

$${\chi }_{\mathrm{sed}}^2=\sum _{i=1}^{N_{\mathrm{sed}}}\frac{{\left(m_{Vi}^{\prime }-m_{Vi}\right)}^2}{{\sigma }_{\mathrm{sed}i}^2}.$$

$${\chi }_{\mathrm{mass}}^2=\sum _{j=1}^{N_{\mathrm{bod}}}{\left(\frac{2m_j-m_j^{\min }-m_j^{\max }}{m_j^{\max }-m_j^{\min }}\right)}^{100}$$

Full text

celestial mechanics — methods: numerical — binaries (including multiple): close, eclipsing — techniques: radial velocities, photometric, interferometric

1 Introduction

Traditional models of eclipsing binaries have to often account for additional external bodies, most importantly as a third light, which makes depths of primary and secondary minima shallower; a light-time effect, causing periodic variations on $O-C$ diagrams; a precession of the argument of periastron $\omega$, shifting the secondary minimum due to perturbations by the 3rd body; or changes of the inclination $i$ with respect to the sky plane, in other words disappearing eclipses.

While analytical theories exist for a description of dynamical perturbations in triple stellar systems and corresponding transit timing variations (also known as TTVs, ETVs; see e.g. Brown 1936, Harrington 1968, Söderhjelm 1975, Breiter & Vokrouhlický 2015, Borkovits et al. 2016), we would prefer a more general approach — to account for all observational data; or at least as many as feasible. So, our aim is to incorporate astrometric or speckle-interferometric positions, radial velocities, minima timings, eclipse durations, spectro-interferometric visibilities, closure phases, synthetic spectra, spectral-energy distribution, and light curves too. At the same time, we do not want to be limited by inevitable approximations of the analytical theories (the N-body problem is not integrable) and the only way out seems to be an N-body integrator (as in Carter et al. 2011).

Another aspect is we cannot use analytical photometric models (like those used for exoplanet transits; Mandel & Agol 2002, Carter et al. 2008, Pál 2012), because the respective simplifications are not acceptable for stellar eclipses, not speaking about ellipsoidal variations outside eclipses.

In principle, our approach should be rather straightforward: we merge two codes into a single one; namely Levison & Duncan (1974) SWIFT code, and Wilson & Devinney (1971) WD code. In practice, a lot of work has to be done, because both of them have to be modified, we need to extract and derive observable quantities, read observational data and check them by means of a $\chi^{2}$ statistics. Last but not least, we need to run a minimisation algorithm on top of them.

Even though we do not present new observational data here, there is one recent application of our N-body model to $\xi$ Tauri quadruple system which was described in a great amount of detail in Nemravová et al. (2016). Moreover, there is a comparison with a number of traditional, observation-specific models. In this ‘technical’ paper, we prefer to show mostly results of numerical simulations, or even negative results contradicting the observations, to demonstrate a sensitivity of our model.

We have a few motivations to do so: (i) no complete and fully self-consistent N-body model exists yet, which can account for that many observational constraints, (ii) we improved the model significantly compared to Nemravová et al. as we can now fit also complete light curves and optionally individual spectra (to be matched by synthetic ones); (iii) the previous paper was a bit lengthy and there was simply not enough room for a more technical description of our code; (iv) we have to discuss the role of systematics, an experience gained during modelling of real multiple stellar systems.

2 Model description

Let us begin with a description of the numerical integrator and the photometric model; then we present principal equations, a definition of the $\chi^{2}$ metric used to compare the model with observational data, and a list of dynamical effects that can be modelled.

2.1 Numerical integrator

We use the Bulirsch–Stoer numerical integrator (Press et al. 1999), with an adaptive time step, controlled by a unit-less parameter $\epsilon_{\rm BS}$. The integrator sequentially divides the time step $\Delta t$ by factors 2, 4, 6, …, checks if the relative difference between successive divisions is less than $\epsilon_{\rm BS}$ and then performs an extrapolation $\Delta t\to 0$ by means of a rational function (see Figure 1). If the maximum number of divisions $n_{\rm max}=10$ is reached, the basic time step $\Delta t$ has to be decreased, with another maximum number of trials $n_{\rm try}=30$. We beg to recall this well-known principle here as it is important to always understand principles and limitations of numerical methods in use. This kind of integrator is quite general and there are no restrictions for magnitudes of perturbations, so we can handle keplerian orbits, tiny N-body perturbations or even violent close encounters. Even though it is not symplectic, it does not suffer from an artificial periastron advance. On long time scales, it is worth to check the energy conservation and eventually decrease $\epsilon_{\rm BS}$, perhaps down to $10^{-11}$.

Apart from the internal time step, a user can choose the output time step $\Delta t_{\rm out}$. The time stepping was adapted so that we first prepare a list of ‘times of interest’ (corresponding to all observations) and the integrator outputs coordinates and velocities at exactly these times. Consequently, the need for additional interpolations is eliminated, except for minima timings and eclipse durations, where a linear interpolation from two close neighbouring points separated by the expected duration is used, and optionally for light curves (see below).

2.2 Photometric model

The only restriction for the geometry of the stellar system, is that only bodies 1 and 2 may be components of an eclipsing binary (or an ellipsoidal variable). Nevertheless, there can be any number of additional bodies, which do contribute to the total light, but we do not compute eclipses for them.

For light curve computations, we use the WD 2005 version, in order to produce compatible and comparable results to Phoebe 1.0 (Prša & Zwitter, 2005), but we plan to upgrade in the future. In brief, the WD code accounts for: black-body radiation or the Kurucz atmospheres, bolometric limb darkening, gravity darkening, reflection, an axial rotation, or the Rossiter–McLaughlin effect. This is a relatively complex photometric model (more complex than analytical models of Mandel & Agol 2002, Carter et al. 2008). We use no spots or circumstellar clouds in this version. Usually, the code is called with mode 0 (no constraints on potentials) or 2 (the luminosity $L_{2}$ of the secondary is computed from the temperature $T_{2}$). Note a number of parameters in lc.in input file are useless (e.g. orbital elements, precession and period rates, luminosities, potentials etc.) because they are driven from elsewhere.

To speed up light curve computations, we can use a binning of times $\Delta t_{\rm bin}$ and then linearly interpolate light curve points to the times of observations. For high-cadence data, we can possibly gain a factor of 10 or 100 speed-up this way, but we have to be sure there is no physical process in our model which could change magnitudes on the timescale shorter than $\Delta t_{\rm bin}$.

2.3 Principal equations

Principal equations of our N-body model can be summarized as follows (the notation is described in Table 1) — the equation of motion:111The program, including sources and example input data, is available at http://sirrah.troja.mff.cuni.cz/~mira/xitau/.

[TABLE]

with the lowest-order tidal term (Hut 1981):

[TABLE]

oblateness:

[TABLE]

and parametrized post-newtonian (PPN; Mardling & Lin 2002) terms:

[TABLE]

Apart from trivial sky-plane positions $x_{{\rm b}i}$, $y_{{\rm b}i}$ and radial velocities $v_{z{\rm b}i}$, we can derive a number of dependent quantities, such as mid-eclipse timings (including light-time effects):

[TABLE]

eclipse durations:

[TABLE]

luminosities (alternatively, using a black-body approximation, $\pi B_{\lambda}(T_{\rm eff})$):

[TABLE]

a limb-darkened complex visibility (Hanbury Brown et al. 1974; $\Theta=\pi\theta_{j}\!\sqrt{u^{2}\!+\!v^{2}}$, $\alpha=1-u_{\rm limb}$, $\beta=u_{\rm limb}$):

[TABLE]

with $u_{\rm limb}(\lambda,T_{{\rm eff}\,j},\log g_{j},{\cal Z}_{j})$ interpolated from Van Hamme 1993); a complex triple product:

[TABLE]

the true phase of the eclipsing binary (at a time $t$ modified by the light-time effects):

[TABLE]

its inclination:

[TABLE]

Kopal potential (for the WD code which outputs relative magnitudes $m_{V}^{\prime}$):

[TABLE]

[TABLE]

a normalized synthetic spectrum (with appropriate Doppler shifts):

[TABLE]

or a spectral-energy distribution (in any of the UBVRIJHK bands):

[TABLE]

[TABLE]

where the component spectra (both $I_{\rm syn}$ and $F_{\rm syn}$) can be either user-supplied, or interpolated on the fly with Pyterpol (Nemravová et al. 2016) from AMBRE, POLLUX, BSTAR, OSTAR or PHOENIX grids (Palacios et al. 2010, de Laverny et al. 2012, Lanz & Hubený 2007, Lanz & Hubený 2003, Husser et al. 2013).

Internally, we use a barycentric left-handed Cartesian coordinate system with $x$ negative in the right-ascension direction, $y$ positive in declination, and $z$ positive in radial, i.e. away from the observer; the units are day, au, au/day and ${\rm au}^{3}/{\rm day}^{2}$ for the time, coordinates, velocities and masses, respectively. We also need additional coordinate systems, namely: Jacobian (for computations of hierarchical orbital elements), 1-centric (for an eclipse detection), 1+2 photocentric, or 1+2+3 photocentric (for a comparison with astrometric observations of components 3 and 4).

One may immediately note a minor caveat of our model: the geometric radius (in Eq. (6)), the effective radius (in Eq. (7)), the limb-darkened radius (a.k.a. $\theta_{j}$ in Eq. (8)), and the average radius (used in Eq. (12)) are all assumed to be approximately the same. If this does not hold, it would be necessary to add some three more equations describing relations between them.

2.4 Observational data

When we compare our model with observations, we can compute $\chi^{2}$ for astrometric positions, radial velocities, minima timings (TTVs), eclipse durations, interferometric squared visibilities, closure phases, triple product amplitudes, light curves, synthetic spectra, and spectral-energy distribution:

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Again, the quantities are described in Table 1. The index $i$ always corresponds to observational data, $j$ to individual bodies, and $k$ to sets of data. The primed quantities correspond to synthetic data, integrated (or interpolated) to the times of observations $t_{i}$.

We can also add an artificial term:

[TABLE]

to keep the masses $m_{j}$ of the components within resonable intervals (e.g. according to spectroscopic classifications of the components). The high exponent of the arbitrary function prevents simplex to drift away from the interval $(m_{j}^{\rm min},m_{j}^{\rm max})$.

As usually, observational data have to be in a suitable format and we provide some example scripts for a conversion or extraction of data from OIFITS files (Pauls et al. 2005). Note that one shall not use RV measurements when it is possible to fit the observed spectra with synthetic ones. Similarly, no minima timings or durations are needed when we have complete light curves at disposal (cf. Figure 2); and no triple product amplitudes $|T_{3}|$ when the same interferometric measurements are used as squared visibilities $|V|^{2}$. Let us emphasize that it is always better to use directly observable quantities, not derived!

To find a local or a global minimum, we can use a standard simplex algorithm or simulated annealing (Nelder & Mead 1965, Press et al. 1999), with the cooling schedule ${\cal T}^{i+1}=(1-\epsilon_{\rm temp}){\cal T}^{i}$, after given number of iterations at ${\cal T}^{i}$. Free parameters of the model (which can be optionally fixed) are: the masses $m_{j}$ of the components, orbital elements $a_{j},e_{j},i_{j},\Omega_{j},\omega_{j},M_{j}$ of the respective orbits, systemic velocity $\gamma$, distance $d$, radii $R_{j}$, effective temperatures $T_{{\rm eff}\,j}$, projected rotational velocities $v_{{\rm rot}\,j}$, and magnitude zero points $m_{0\,k}$. For $N_{\rm bod}$ bodies, this represents a set of $({\bf 10}N_{\rm bod}+N_{\rm band}-4)$ parameters in total.

Unfortunately, neither of the numerical methods can guarantee that the global minimum will be found. The multidimensional parameter space is very extended and there are many local minima, some of them statistically equivalent. Running simplex many times from different starting points ($10^{4}$ to $10^{5}$) can help, but this problem is clearly both system-dependent and data-dependent. To obtain uncertainties of the model parameters or full covariance matrix one can use the bootstrap method (Efron 1979) for example.

2.5 Dynamical effects

A number of well-known dynamical effects can be modelled with the N-body integrator: self-consistent precession of $\omega$ and $\Omega$, inclination changes and eclipse durations, eccentricity oscillations (Figure 3), Kozai cycles, variation and evection, differences between prograde vs retrograde orbits, close encounters, hyperbolic trajectories, mean-motion resonances (Rivera et al. 2005), secular resonances, three-body resonances (Nesvorný and Morbidelli 1998), or chaotic diffusion due to overlapping resonances are also naturally accounted for in our N-body model. Even more examples can be found in Fabrycky (2010).

3 Model testing

3.1 A test on synthetic data

In order to test a basic functionality of our N-body model, we created a mock system with 40 known parameters; the system actually closely corresponds to the quadruple star $\xi$ Tau (see Table 2). Synthetic observational data were created using the same code with the same coverage and cadence as the real observations of $\xi$ Tau (Nemravová et al. 2016): 78 133 spectral measurements (individual data points $I_{\lambda\,i}$), 17 391 squared visibilities $|V|^{2}_{i}$, 4 856 complex triple products $T_{3\,i}$, 2 974 lightcurve points $m_{Vki}$, 17 astrometric measurements $\Delta x_{ij}$, $\Delta y_{ij}$ of the 4th component, and 13 SED points $m_{Vi}$. A gaussian noise was applied to all of them, at the levels usual for datasets we have for $\xi$ Tau; and we assumed there are no systematics, neither in these synthetic observations, nor in the model (but cf. Appendices A.3 to A.8). For the true solution, one would obtain $\chi^{2}=109\,095$, which is indeed a perfect solution given the number of degrees of freedom $\nu\equiv N_{\rm data}-M_{\rm free}=108\,257-40=108\,217$ and the probability $P(\chi^{2}|\nu)=0.970$ that the $\chi^{2}$ value is that large by chance.

We then performed several simplex optimisations, starting from a number of neighbouring points, farther and farther away from the true solution. The convergence of the simplex to a local minimum is clear (Figure 4) and we recovered the original parameters (low $\chi^{2}\simeq\nu$) with some uncertainties, as expected; unless the initial guess was too far away, say more than a few percent of the critical parameters (see Figure 5). More extended surveys with many initial starting points and/or simulated annealing would be needed in these unfortunate cases. The differences between the final and true solutions roughly correspond to the uncertainties we would obtain from bootstrap testing.

Often some preliminary knowledge based on observation-specific models is available (e.g. prominent periods, previously published parameters). The N-body model is especially suitable for such ‘semi-final’ convergence — with all parameters free — which should limit the systematics arising from a usage of limited (keplerian) models.

Regarding the fractional or missing data, there are several rather trivial facts, e.g. if we miss RVs, it is impossible to resolve low-$e$ orbits from high-$e$ with $\omega=0^{\circ}$ or $180^{\circ}$. If there are no eclipses and no interferometric measurements available, one cannot precisely constrain the inclination $i$. Without closure phase measurements, there is practically no sensitivity to asymmetries, etc. Of course, they are interesting when dealing with real observational datasets (see also Appendix A.2).

Let us point out that this kind of testing has somewhat limited capabilities. First, if we create the mock data with the same model, then this test is essentially the test of the numerical methods (simplex and simulated annealing). Because these methods are indeed classical (Nelder & Mead 1969), their limitations are already very well known.

Second, to create the mock data independently, one would need a completely independent model with exactly the same capabilities. However, this is rather a test of systematic differences between the models and not of the method itself. We consider this approach to be more useful, but it is rather difficult to obtain the second model. Of course, any keplerian models (e.g. Phoebe 1.0) are useless. Any model which does not produce all the observables (astrometry, RVs, TTVs, $\epsilon_{\rm ecl}$, $|V|^{2}$, $T_{3}$, $m_{V}$, or $I_{\lambda}$) is impractical too, because we need as much orthogonal constraints as possible. Analytical photometric models (e.g. Carter et al. 2008) are too simplified for stellar binaries. And so on. It may be possible to use Phoebe 2.0 (Prša et al. 2016) for such comparison in the future.

Last but not least, real observational data of real systems have its own cadence, coverage, calibrations, uncertainties, and systematics. Even though one can play with artificial data, these tests cannot be used straightforwardly, because we know ”nothing” a-priori about the given data we obtain from observers. Consequently, one will have to perform suitable tests again and again for every next system.

3.2 A comparison with other models

As already mentioned in the Introduction, a comparison with several observation-specific models was already done in our previous paper Nemravová et al. (2016). In particular, our N-body model produces results which are compatible within respective approximations (e.g. on short time scales when the keplerian model can be regarded as a useful approximation) with the following published works: photometric model of Phoebe 1.0 (Prša & Zwitter 2006); astrometric and speckle-interferometry model of Zasche & Wolf (2007); RV disentangling by Korel (Hadrava 1997); LitPro model for visibilities $|V|^{2}$ (Tallon-Bosc et al. 2008); spectro-interferometric $|V|^{2}$ and $T_{3}$ (Nemravová et al. 2016); synthetic spectra fitting by Pyterpol (dtto).

4 Conclusions and future work

Today, N-body models seem to be an absolutely necessary tool for a careful inspection of observational data. It is important to take care that discrepancies between keplerian and full N-body dynamics no longer spoil derived stellar parameters. After a removal of (some) systematic errors (sometimes) present in observations or reductions, it enables us to reveal even tiny N-body perturbations and construct robust models of compact stellar systems (e.g. those from Table 3).

Regarding future developments of (our or other) N-body models, it seems worthwhile to also account for: calibration factors of individual interferometric telescopes, gravity darkening in the visibility calculation of rotating stars (as in Aufdenberg et al. 2006), especially when measuring on longest baselines, and eventually one may think of an upgrade to the WD 2015, or Phoebe 2.0 (already used in Pablo et al. 2015).

Yet another work is needed to compute trajectories even more accurately, with physics going beyond point-like masses, equilibrium tides or oblateness, namely: higher gravitational moments ($J_{4}$) due to non-sphericity of stellar components, tidal dissipation and cross tides (e.g. Mignard 1979), corresponding long-term evolution of orbits, spin evolution (Eggleton & Kisileva-Eggleton 2001), spin–orbital resonances, or radiation of gravitational waves in extreme cases.

The situation in stellar interiors also matters. The dissipation occurs either due to viscosity in outer convective zones, or due to inertial oscillations in radiative zones, which are excited on eccentric orbits by dynamic tides and subsequently radiatively damped (Zahn 2008). In triple systems, the excitations may actually arise from a binary subsystem, and corresponding light oscillations then have an half of its period (Derekas et al. 2011, Fuller et al. 2013). Another difficulty stems from certain coupling of envelopes and cores (Papaloizou & Ivanov 2010). Inevitably, a fully self-consistent model should account for a back-reaction: the strongest tidal heating may inflate whole objects (Mardling 2007).

The work of MB was supported by the grants no. P209-15-02112S and P209-13-01308S of the Czech Science Foundation. I thank Jana Nemravová and David Vokrouhlický for valuable discussions on the subject and a fruitful collaboration on the $\xi$ Tauri paper. I also have to thank the referee Hagai Perets for a constructive criticism which determined the final structure of the paper.

Appendix A Possible problems due to systematics

We have to admit that any modelling (compact stellar systems included) can be spoiled, either when there are systematic deficiencies of the model, e.g. keplerian vs N-body, or serious systematic errors in observational data, especially when we use very heterogeneous datasets. In the following, we thus discuss several ‘dangerous’ cases.

A.1 Discretization errors

Of course, any numerical computation suffers from discretization errors and interpolation errors, even though we tried to decrease the latter as much as possible (cf. Section 2). This is probably the most important disadvantage compared to analytical computations. A general rule is a convergence of results (and corresponding $\chi^{2}$ values) for $\Delta t\to 0$.

However, let us add a warning that rarely a decrease of time step, e.g. by a factor of 2, may lead to unexpected results. For example, when eclipses are almost disappearing, the trajectory with $\Delta t/2$ is more curved and may thus miss the last eclipse, which suddenly increases $\chi^{2}_{\rm ttv}$ because the next eclipse is now one orbital period $P$ far away. The solution is to converge the model once again, with $\Delta t/2$.

Let’s not forget, there is yet another discretization related to the WD code, or the surfaces of the eclipsing binary. For low numbers $N_{\rm wd}$, one can see numerical artefacts on the light curve, as rectangular surface facets appear from behind the limb, or disappear. Again, it is worth to check larger $N_{\rm wd}$.

A.2 Mirror solutions

Quite often, we can expect one or more ($m$) mirror solutions (and $2^{m}$ combinations of them). A typical situation is we have no RVs for faint components (so that both inclinations $i$ and $i^{\prime}=-i$ are admissible), or no unambiguous astrometry or closure-phase measurements (so that $\Omega$ and $\Omega^{\prime}=180^{\circ}-\Omega$ are both admissible). Consequently, one may save some time when surveying the parameter space.

However, with the N-body model at hand it is worth to check not only the total $\chi^{2}$ but also individual contributions to $\chi^{2}$ for all the mirror models! Especially $\chi^{2}_{\rm ttv}$ is very sensitive to the mutual perturbations, and we may be able to resolve some of the ambiguities mentioned above.

Of course, the statistics must not be corrupted by systematics or strongly underestimated uncertainties in other observational datasets. If this is the unfortunate case, one may try to use weights $w$ of individual $\chi^{2}$’s, but this should be used as “a method of last resort”. The reason is that it is too easy to hide all systematics this way, even though it is better to get rid of them (see below).

A.3 Heterogeneous datasets of RVs

Radial-velocity measurements might be affected by zero point offsets, which then lead to different systemic velocities $\gamma$ for different observatories. This can be a bit misleading, because it is not possible to a priori distinguish systematic differences in dispersion relations from real perturbations, when the observations were acquired at epochs distant in time.

A well-known viable approach is to use an independent calibration by narrow interstellar lines (DIBs; Chini et al. 2012), if they are present and resolved in the given spectral range. Another possibility are atmospheric lines for which the relative RVs can be computed easily. If this is impossible, one should use the N-body model with a great caution, because simply increasing $\sigma_{\rm rv}\simeq\Delta\gamma$, to get $\chi^{2}_{\rm rv}\simeq N_{\rm rv}$ is a wrong idea. The RV measurements in question will still ‘push’ the model elsewhere and there will be systematic departures with respect to other (more or less orthogonal) observational data.

It may be a too much freedom, but if the dispersion relations can be considered stable from night to night, some calibration factors $f_{{\rm rv}\,k}$ of RVs — assigned to individual observatories or datasets — might be actually a better solution. In any case, such factors have to be always treaded as additional free parameters of the N-body model.

A.4 RVs from disentangling

Sometimes, RVs are derived in the Fourier domain by means of disentangling (e.g. by Korel; Hadrava 1995), with an advantage to obtain disentangled spectra of individual components. There is a ‘hidden’ caveat, though, because one can expect a strong correlation of RVs and the fixed keplerian orbital elements used during the disentangling procedure. This represents a problem, because we do vary initial osculating orbital elements in the N-body model and they most likely will contradict the previous elements.

Note the disentangled spectra should not be re-used as templates, because they contain slight systematic asymmetries or wavy continua. If we try to match the observed spectra with such templates again, we would obtain artificially small uncertainties $\sigma_{\rm rv}$ (and extremely large $\chi^{2}_{\rm rv}$). A solution is to use synthetic spectra similar to the disentangled ones, but with no direct relation to Korel, as an intermediate step to derive new RVs.

A.5 RVs from synthetic spectra

Alternatively, RVs of the individual components can be derived directly in the time domain by fitting a luminosity-weighted sum of suitable synthetic spectra (e.g. by Pyterpol; Nemravová et al. 2016). Instead of fitting the observed spectra individually (one-by-one), it is advisable to assume that most of the free parameters (projected $v_{{\rm rot}\,j}$, $T_{{\rm eff}\,j}$, gravity $\log g_{j}$, and metallicity $Z_{j}$ of the stellar components) are the same for all spectra, with the exception of RVs which are surely time dependent. Luckily, these RVs are not strongly correlated with the orbital elements, so they seem suitable as an input for the N-body model.

On the other hand, this method can have problems on its own when RVs are small (at conjunctions) and $v_{\rm rot}$ large, so that the lines are totally blended. As a provisional solution, one may try to discard the lowest RVs which cause the problems, or do not use RVs at all and rather fit synthetic spectra directly with the N-body model ($\chi^{2}_{\rm syn}$), which is definitely a better approach, because RVs will be correctly tied to each other (see Figure 6).

A.6 Rectification procedure

Inevitably, RVs might have been systematically affected already during a basic reduction, namely a rectification (normalisation) of spectra. If the rectification procedure is automated by fitting a low-degree polynomial to continua, it is worth to try a different maximum degree of the polynomial and run the above synthetic spectra optimisation once again.

A.7 Visibility calibration

Contrary to closure phase $\arg T_{3}$ measurements, the squared visibility $|V|^{2}$ has to be calibrated by close-in-time observations of comparison stars with known angular diameters or unresolved (point-like) sources. Sometimes even the calibrated measurements exhibit unrealistically quick changes of $|V|^{2}$ or sudden decreases of $|V|^{2}$, possibly caused by unfavorable weather conditions, or seeing comparable to the slit width, affecting a light contribution from barely-resolved components, or other obscure instrumental defects.

In the end, dropping of these suspicious observational data may be the only way to prevent the systematics to unrealistically shift the model. Using a low weight $w_{\rm vis}=0.1$ is not a satisfactory option. To this point, we always retain a dataset identification for each single measurement which enables us to quickly perform a bootstrap testing.

A.8 Quasiperiodic oscillations

A removal of quasiperiodic light oscillations which are sometimes present (or always for high-precision measurements) outside eclipses is very important, because they may otherwise systematically offset the minima timings themselves. One wave of the oscillations behaves like a ‘ramp’, which skews the light curve at around the minimum.

The observed light curve should be thus locally fitted by a suitable function (e.g. harmonic with a variable period and amplitude) and then subtracted from the data. If the (synthetic) light curve out of eclipses is flat beyond doubt, it seems better to drop these segments of the (observed) light curve completely, because they would increase $\chi^{2}_{\rm lc}$ but there is no useful information as we have no physical model for these oscillations (as of yet).

A.9 Osculating vs fixed elements

Some care is also needed when comparing results of (old) keplerian and (new) N-body models. They actually can differ by more than a few $\sigma$, because the former orbital elements are fixed, while the latter are only osculating initial conditions at $t=T_{0}$. Generally, all elements are time-dependent quantities, $a_{1}(t),e_{1}(t),i_{1}(t),$ etc., whereas their oscillations are often larger than uncertainties of the initial osculating elements. In fact, one can perform some averaging over the observational time span to facilitate the comparison. Nevertheless, the N-body model is more complete, and it should be probably preferred.

A.10 Stability, aliasing, mean and proper elements

It is also possible to run the N-body integrator separately, regardless of an observational time span, and study a long-term evolution and stability of stellar systems. We may wish to prefer those orbital solutions which are indeed stable. One of the difficulties is that the output of osculating elements is either prohibitively long or an aliasing occurs when the output time step $\Delta t_{\rm out}$ is larger than an half of the shortest orbital period, $P_{1}/2$.

In a modified version of the BS integrator (swift_bs_fp), we can use an on-line digital filtering of non-singular osculating elements $h_{j},k_{j},p_{j},q_{j}$ to overcome these problems: first a multi-level convolution based on the Kaiser windows (Quinn et al. 1991) to obtain mean elements, and second a frequency-modified Fourier transform (Šidlichovký & Nesvorný 1997) to extract proper elements. For $N$ mutually interacting bodies, one can expect $2N$ eigen-frequencies of the system, which are usually denoted $g_{j}$ and $s_{j}$. The corresponding amplitudes $e_{{\rm p}j}$, $\sin{1\over 2}I_{{\rm p}j}$ can be considered approximate integrals of motion which only evolve on time scales longer than secular.

To conclude in a pessimistic way, the above list of possible problems and systematics cannot be treated as complete, unfortunately.

Appendix B Technical notes

B.1 Different hierarchy

By default, we assume a hierarchy like ((1+2)+3)+4, for which Jacobian orbital elements seem to be a suitable description. For a substantially different hierarchy, say two pairs like (1+2) and (3+4), where we would prefer a different definition of elements, only a very small part of the code has to be rewritten, namely in the geometry.f subroutine, where the elements are converted to barycentric Cartesian coordinates. Alternatively, one may wish to use 1-centric Cartesian coordinates as actual parameters, because sometimes precise observational data may constrain (‘fix’) some of them ($x_{{\rm h}2}$, $y_{{\rm h}2}$, etc.), thus decreasing the dimensionality of the parameter space.

B.2 Jacobian orbital elements

Unlike the usual stellar-astronomy convention, where the brightest component is always at the origin of the reference frame, in our N-body model we usually select the most compact eclipsing pair as bodies 1 and 2, or the most massive component as 1. The reason is that orbital elements in hierarchical systems are usually computed in Jacobian coordinates, where the centre of mass 1+2 is the reference point for the coordinates and velocities of the 3rd body; the 1+2+3 centre of mass is a suitable reference for the 4th body, and so on. The corresponding Jacobian elements then have a nice interpretation. Because of the above definition, it may be necessary to adjust to-be-fitted astrometric measurements by $180^{\circ}$ in the position angle — not due to an ambiguity, but simply because the reference body is different in our case. Similarly, a value of $\Omega$ from literature may actually differ by $180^{\circ}$.