An Advanced N-body Model for Interacting Multiple Stellar Systems

We use the Bulirsch–Stoer numerical integrator (Press et al. 1999), with an adaptive time step, controlled by a unit-less parameter _BS. The integrator sequentially divides the time step Δt by factors 2, 4, 6, ..., checks if the relative difference between successive divisions is less than _BS and then performs an extrapolation Δt → 0 by means of a rational function (see Figure 1). If the maximum number of divisions n_max = 10 is reached, the basic time step Δt has to be decreased, with another maximum number of trials n_try = 30. We recall this well-known principle here, as it is important to always understand the principles and limitations of the 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 timescales, it is worthwhile to check the energy conservation and eventually decrease _BS, perhaps down to 10⁻¹¹.

Zoom In Zoom Out Reset image size

Apart from the internal time step, a user can choose the output time step Δt_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 neighboring points separated by the expected duration is used, and optionally one can use binning for light curves (see below).

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: blackbody 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₂ of the secondary is computed from the temperature T₂). Note that a number of parameters in the 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 Δt_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 that could change magnitudes on a timescale shorter than Δt_bin.

Principal equations of our N-body model can be summarized as follows (the notation is described in Table 1): the equation of motion¹

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{r}}}}_{{\rm{b}}i}=-\displaystyle \sum _{j\ne i}^{{N}_{\mathrm{bod}}}\displaystyle \frac{{{Gm}}_{j}}{{r}_{{ji}}^{3}}\,{{\boldsymbol{r}}}_{{ji}}+{{\boldsymbol{f}}}_{\mathrm{tidal}}+{{\boldsymbol{f}}}_{\mathrm{oblat}}+{{\boldsymbol{f}}}_{\mathrm{ppn}}\,;\end{eqnarray} \tag{ 1 }$$

with the lowest-order tidal term (Hut 1981)

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{tidal}}=-\displaystyle \sum _{j\ne i}^{{N}_{\mathrm{bod}}}3{k}_{{\rm{L}}i}\displaystyle \frac{{{Gm}}_{j}^{2}}{{m}_{i}}\displaystyle \frac{{R}_{i}^{5}}{{r}_{{ij}}^{8}}{{\boldsymbol{r}}}_{{ji}};\end{eqnarray} \tag{ 2 }$$

oblateness

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{oblat}}=-\displaystyle \sum _{j\ne i}^{{N}_{\mathrm{bod}}}\displaystyle \frac{1}{2}{k}_{{\rm{L}}j}\,{\omega }_{\mathrm{rot}\,j}^{2}\displaystyle \frac{{R}_{i}^{5}}{{r}_{{ij}}^{5}}{{\boldsymbol{r}}}_{{ji}};\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}{{\boldsymbol{f}}}_{\mathrm{ppn}} & = & -\,\displaystyle \sum _{j\gt i}^{{N}_{\mathrm{bod}}}\displaystyle \frac{G({m}_{j}+{m}_{i})}{{r}_{{ij}}^{2}{c}^{2}}\phantom{\rule{}{3.5ex}}\{-2(2-{\eta }_{{ij}}){\dot{r}}_{{ij}}{\dot{{\boldsymbol{r}}}}_{{ij}}\\ & & +\left[(1+3{\eta }_{{ij}}){\dot{{\boldsymbol{r}}}}_{{ij}}\cdot {\dot{{\boldsymbol{r}}}}_{{ij}}-\displaystyle \frac{3}{2}{\eta }_{{ij}}{\dot{r}}_{{ij}}^{\,2}\right.\\ & & -\left.\left.2(2+{\eta }_{{ij}})\displaystyle \frac{G({m}_{j}+{m}_{i})}{{r}_{{ij}}}\right]\displaystyle \frac{{{\boldsymbol{r}}}_{{ij}}}{{r}_{{ij}}}\right\}.\end{eqnarray} \tag{ 4 }$$

Apart from trivial sky-plane positions x_bi, y_bi and radial velocities v_zbi, we can derive a number of dependent quantities, such as mid-eclipse timings (including light-time effects)

$$\begin{eqnarray}&&{t}_{\mathrm{ecl}}^{\prime }={t}_{\min }+\displaystyle \frac{{z}_{{\rm{b}}1+2}-{z}_{{\rm{b}}1+2}(t={T}_{0})}{c}-\displaystyle \frac{{z}_{{\rm{h}}2}}{c};\end{eqnarray} \tag{ 5 }$$

eclipse durations

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ecl}}^{\prime }=\displaystyle \frac{2}{{\bar{v}}_{{\rm{h}}2}}\sqrt{{({R}_{1}+{R}_{2})}^{2}-{{\rm{\Delta }}}_{\min }^{2}};\end{eqnarray} \tag{ 6 }$$

luminosities (alternatively, using a blackbody approximation, πB_λ(T_eff))

$$\begin{eqnarray}{L}_{j}({T}_{\mathrm{eff}j},{R}_{j}) & = & 4\pi {R}_{j}^{2}\\ & & \times {\displaystyle \int }_{\lambda -{\rm{\Delta }}\lambda /2}^{\lambda +{\rm{\Delta }}\lambda /2}{F}_{\mathrm{syn}}\,(\lambda ,{T}_{\mathrm{eff}j},\mathrm{log}{g}_{j},{v}_{\mathrm{rot}j},{{ \mathcal Z }}_{j})d\lambda ;\end{eqnarray} \tag{ 7 }$$

a limb-darkened complex visibility (Hanbury Brown et al. 1974; ${\rm{\Theta }}=\pi {\theta }_{j}\sqrt{{u}^{2}+{v}^{2}}$, α = 1 − u_limb, β = u_limb)

$$\begin{eqnarray}V^{\prime} (u,v) & = & \displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}\displaystyle \frac{{L}_{j}}{{L}_{\mathrm{tot}}}{\left(\displaystyle \frac{\alpha }{2}+\displaystyle \frac{\beta }{3}\right)}^{-1}\\ & & \times \left[\alpha \displaystyle \frac{{J}_{1}({\rm{\Theta }})}{{\rm{\Theta }}}+\,\beta \sqrt{\displaystyle \frac{\pi }{2}}\displaystyle \frac{{J}_{3/2}({\rm{\Theta }})}{{{\rm{\Theta }}}^{3/2}}\right]\,{e}^{-2\pi i({{ux}}_{{\rm{a}}j}+{{vy}}_{{\rm{a}}j})};\end{eqnarray} \tag{ 8 }$$

with ${u}_{\mathrm{limb}}(\lambda ,{T}_{\mathrm{eff}j},\mathrm{log}{g}_{j},{{ \mathcal Z }}_{j})$ interpolated from Van Hamme 1993); a complex triple-product

$$\begin{eqnarray}&&{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}));\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}&&{\varphi }_{\mathrm{ecl}}^{\prime }=\displaystyle \frac{1}{2\pi }\arctan \displaystyle \frac{\hat{O}\cdot \hat{Y}}{\hat{O}\cdot \hat{X}};\end{eqnarray} \tag{ 10 }$$

its inclination

$$\begin{eqnarray}&&{i}_{\mathrm{ecl}}^{\prime }=\arccos (-\hat{O}\cdot \hat{Z});\end{eqnarray} \tag{ 11 }$$

Kopal potential (for the WD code that outputs relative magnitudes ${m}_{V}^{\prime }$)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{Kopal}j}\doteq {\left\langle \displaystyle \frac{1}{{r}_{1}}+\displaystyle \frac{q}{{r}_{2}}+\displaystyle \frac{1}{2}(1+q){r}_{3}^{2}\right\rangle }_{\mathrm{circle}\displaystyle \frac{{R}_{j}}{{r}_{12}},\mathrm{where}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray*}&&{{\boldsymbol{r}}}_{1}\in \mathrm{circle},{{\boldsymbol{r}}}_{2}={{\boldsymbol{r}}}_{1}-(1,0,0),{{\boldsymbol{r}}}_{3}=\left({x}_{1}-\displaystyle \frac{q}{1+q},{y}_{1},0\right);\end{eqnarray*}$$

a normalized synthetic spectrum (with appropriate Doppler shifts)

$$\begin{eqnarray}&&{I}_{\lambda }^{\prime }=\displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}\displaystyle \frac{{L}_{j}}{{L}_{\mathrm{tot}}}\,{I}_{\mathrm{syn}}\,\left[\lambda \left(1-\displaystyle \frac{{v}_{z{\rm{b}}j+\gamma }}{c}\right),{T}_{\mathrm{eff}j},\mathrm{log}{g}_{j},{v}_{\mathrm{rot}j},{{ \mathcal Z }}_{j}\right];\end{eqnarray} \tag{ 13 }$$

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

$$\begin{eqnarray}&&{F}_{V}^{\prime }=\displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}{\left(\displaystyle \frac{{R}_{j}}{d}\right)}^{2}{\int }_{0}^{\infty }{F}_{\mathrm{syn}}\,[\lambda ,{T}_{\mathrm{eff}j},\mathrm{log}{g}_{j},{v}_{\mathrm{rot}j},{{ \mathcal Z }}_{j}]\,{f}_{V}(\lambda )d\lambda ,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray*}&&{m}_{V}^{\prime }=-2.5{\mathrm{log}}_{10}\displaystyle \frac{{F}_{V}^{\prime }}{{F}_{V\mathrm{calib}}{\int }_{0}^{\infty }{f}_{V}(\lambda )d\lambda },\end{eqnarray*}$$

where the component spectra (both I_syn and F_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 (Lanz & Hubený 2003, 2007; Palacios et al. 2010; de Laverny et al. 2012; Husser et al. 2013).

Table 1.  Notation Used for Coordinates, Velocities, and a Number of Other Quantities and Uncertainties That We Use in Our N-body Model

Nbod	number of bodies
m	mass (GM⊙ units)
$q=\tfrac{{m}_{1}}{{m}_{2}}$	mass ratio
${\eta }_{{ij}}=\tfrac{{m}_{j}{m}_{i}}{{({m}_{j}+{m}_{i})}^{2}}$	symmetrized mass ratio
kL	Love number
ωrot	rotational angular velocity
xb, yb, zb	barycentric coordinates
vxb, vyb, vzb	barycentric velocities
xh, yh, zh	1-centric coordinates
vxh, vyh, vzh	1-centric velocities
xp, yp	1 + 2 photocentric sky-plane coordinates
xp3, yp3	1 + 2 + 3 photocentric coordinates
${x}_{{\rm{a}}}=\tfrac{{x}_{{\rm{h}}}}{d},{y}_{{\rm{a}}}$	1-centric coordinates in an angular measure
$\hat{X},\hat{Y},\hat{Z}$	unit vectors aligned with 1 + 2 eclipsing pair
$\hat{O}=(0,0,-1)$	observers direction
γ	systemic velocity
vrad	observed radial velocity
tecl	mid-epoch of an eclipse of 1 + 2 pair
ecl	eclipse duration
L, Ltot	component luminosity and the total one
Teff	effective temperature
R	stellar radius
λ, Δλ	effective wavelength and bandwidth
Bλ(T)	the Planck function
V	complex visibility; squared visibility is $| V{| }^{2}$
T3	complex triple-product; closure phase is argT3
u, v	projected baselines (expressed in cycles, $\tfrac{B}{\lambda }$)
$\theta =\tfrac{2R}{d}$	angular diameter
ulimb	linear limb-darkening coefficient
d	distance to the system
mV	magnitude (in V band or another)
m0	zero-point
Iλ, Isyn	normalized monochromatic intensity
Fsyn	absolute monochromatic flux (in erg s−1 cm−2 cm−1)
FVcalib	calibration flux
fV	filter transmission coefficient
$g=\tfrac{{GM}}{{R}^{2}}$	surface gravity, $\mathrm{log}g$ in cgs
vrot	projected rotational velocity
${ \mathcal Z }$	metallicity
${\sigma }_{\mathrm{sky}\mathrm{major},\mathrm{minor}}$	uncertainty of the astrometric position,
	angular sizes of the uncertainty ellipse
ϕellipse	position angle of the ellipse
R(...)	the corresponding 2 × 2 rotation matrix
σrv	uncertainty of the radial velocity
σttv	uncertainty of the eclipse mid-epoch timing
σecl	uncertainty of the eclipse duration
σvis	uncertainty of the squared visibility
σclo	uncertainty of the closure phase
σt3	uncertainty of the triple-product amplitude
σlc	uncertainty of the light-curve data
σsyn	uncertainty of the normalized intensity
σsed	uncertainty of the spectral energy distribution
${m}_{j}^{\min },{m}_{j}^{\max }$	minimum and maximum masses

Download table as:  ASCIITypeset image

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 au³/day² 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 Equation (6)), the effective radius (in Equation (7)), the limb-darkened radius (i.e., θ_j in Equation (8)), and the average radius (used in Equation (12)) are all assumed to be approximately the same. If this does not hold, it would be necessary to add three more equations describing the relations between them.

When we compare our model with observations, we can compute χ² 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:

$$\begin{eqnarray}{\chi }^{2} & = & \underline{{\chi }_{\mathrm{sky}}^{2}}+{\chi }_{\mathrm{rv}}^{2}+{\chi }_{\mathrm{ttv}}^{2}+{\chi }_{\mathrm{ecl}}^{2}+\underline{{\chi }_{\mathrm{vis}}^{2}}+\underline{{\chi }_{\mathrm{clo}}^{2}}+{\chi }_{{\rm{t}}3}^{2}+\\ & & +\,\underline{{\chi }_{\mathrm{lc}}^{2}}+\underline{{\chi }_{\mathrm{syn}}^{2}}+\underline{{\chi }_{\mathrm{sed}}^{2}}+{\chi }_{\mathrm{mass}}^{2},\end{eqnarray} \tag{ 15 }$$

where:

$$\begin{eqnarray}&&{\chi }_{\mathrm{sky}}^{2}=\displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{sky}j}}\left\{\displaystyle \frac{{({\rm{\Delta }}{x}_{{ji}})}^{2}}{{\sigma }_{\mathrm{sky}\,\mathrm{major}\,{ji}}^{2}}+\displaystyle \frac{{({\rm{\Delta }}{y}_{{ji}})}^{2}}{{\sigma }_{\mathrm{sky}\,\mathrm{minor}\ {ji}}^{2}}\right\},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&({\rm{\Delta }}{x}_{{ji}},{\rm{\Delta }}{y}_{{ji}})={\boldsymbol{R}}\left(-{\phi }_{\mathrm{ellipse}}-\displaystyle \frac{\pi }{2}\right)\times \left(\begin{array}{c}{x}_{{\rm{p}}\,{ji}}^{\prime }-{x}_{{\rm{p}}{ji}}\\ {y}_{{\rm{p}}\,{ji}}^{\prime }-{y}_{{\rm{p}}{ji}}\end{array}\right),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{rv}}^{2}=\displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{rv}j}}\displaystyle \frac{{({v}_{z{\rm{b}}{ji}}^{\prime }+\gamma -{v}_{\mathrm{rad}{ji}})}^{2}}{{\sigma }_{\mathrm{rv}\,{ji}}^{2}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{ttv}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{ttv}}}\displaystyle \frac{{({t}_{\mathrm{ecl}i}^{\prime }-{t}_{\mathrm{ecl}i})}^{2}}{{\sigma }_{\mathrm{ttv}\,i}^{2}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{ecl}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{ecl}}}\displaystyle \frac{{({\epsilon }_{\mathrm{ecl}i}^{\prime }-{\epsilon }_{\mathrm{ecl}i})}^{2}}{{\sigma }_{\mathrm{ecl}\,i}^{2}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{vis}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{vis}}}\displaystyle \frac{{({| V^{\prime} ({u}_{i},{v}_{i})| }^{2}-| V{| }_{i}^{2})}^{2}}{{\sigma }_{\mathrm{vis}\,i}^{2}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{clo}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{clo}}}\displaystyle \frac{{(\arg {T}_{3i}^{\prime }-\arg {T}_{3i})}^{2}}{{\sigma }_{\mathrm{clo}\,i}^{2}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\chi }_{{\rm{t}}3}^{2}=\displaystyle \sum _{i=1}^{{N}_{{\rm{t}}3}}\displaystyle \frac{{(| {T}_{3i}| ^{\prime} -| {T}_{3i}| )}^{2}}{{\sigma }_{{\rm{t}}3\,i}^{2}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{lc}}^{2}=\displaystyle \sum _{k=1}^{{N}_{\mathrm{band}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{lc}k}}\displaystyle \frac{{({m}_{{Vki}}^{\prime }+{m}_{0k}-{m}_{{Vki}})}^{2}}{{\sigma }_{\mathrm{lc}\,{ki}}^{2}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{syn}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{syn}}}\displaystyle \frac{{({I}_{\lambda i}^{\prime }-{I}_{\lambda i})}^{2}}{{\sigma }_{\mathrm{syn}\,i}^{2}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\chi }_{\mathrm{sed}}^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{sed}}}\displaystyle \frac{{({m}_{{Vi}}^{\prime }-{m}_{{Vi}})}^{2}}{{\sigma }_{\mathrm{sed}\,i}^{2}}.\end{eqnarray} \tag{ 26 }$$

Again, the quantities are described in Table 1. The index i always corresponds to observational data, j corresponds to individual bodies, and k corresponds 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,

$$\begin{eqnarray}&&{\chi }_{\mathrm{mass}}^{2}=\displaystyle \sum _{j=1}^{{N}_{\mathrm{bod}}}{\left(\displaystyle \frac{2{m}_{j}-{m}_{j}^{\min }-{m}_{j}^{\max }}{{m}_{j}^{\max }-{m}_{j}^{\min }}\right)}^{100},\end{eqnarray} \tag{ 27 }$$

to keep the masses m_j of the components within reasonable intervals (e.g., according to spectroscopic classifications of the components). The high exponent of the arbitrary function prevents simplex from drifting away from the interval $({m}_{j}^{\min },{m}_{j}^{\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 should 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 available (cf. Figure 2); and no triple-product amplitudes $| {T}_{3}|$ are needed when the same interferometric measurements are used as squared visibilities $| V{| }^{2}$. We emphasize that it is always better to use directly observable quantities rather than derived ones.

Zoom In Zoom Out Reset image size

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 ${{ \mathcal T }}^{i+1}=(1-{\epsilon }_{\mathrm{temp}}){{ \mathcal T }}^{i}$, after a given number of iterations at ${{ \mathcal 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, Ω_j, ω_j, M_j of the respective orbits, systemic velocity γ, distance d, radii R_j, effective temperatures ${T}_{\mathrm{eff}j}$, projected rotational velocities v_{rot j}, and magnitude zero-points m_{0 k}. For N_bod bodies, this represents a set of $(10{N}_{\mathrm{bod}}+{N}_{\mathrm{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⁴–10⁵) can help, but this problem is clearly both system-dependent and data-dependent. To obtain uncertainties of the model parameters or a full covariance matrix, one can use the bootstrap method (Efron 1979), for example.

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

Zoom In Zoom Out Reset image size

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 ξ Tau (see Table 2). Synthetic observational data were created using the same code with the same coverage and cadence as the real observations of ξ Tau (Nemravová et al. 2016): 78,133 spectral measurements (individual data points I_{λ i}), 17 391 squared visibilities $| V{| }_{i}^{2}$, 4 856 complex triple products ${T}_{3i}$, 2 974 light-curve points m_Vki, 17 astrometric measurements Δx_ij, Δy_ij of the 4th component, and 13 SED points m_Vi. A Gaussian noise was applied to all of them, at levels typical for the data sets we have for ξ Tau; and we assumed there are no systematics, neither in these synthetic observations, nor in the model (but cf. Appendices A.3–A.8). For the true solution, one would obtain χ² = 109,095, which is indeed a perfect solution given the number of degrees of freedom ν ≡ N_data − M_free = 108,257–40 = 108,217, and the probability $P({\chi }^{2}| \nu )=0.970$ that the χ² value is that large by chance.

We then performed several simplex optimizations, starting from a number of neighboring 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 χ² ≃ ν) 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 the 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 ω = 0° or 180°. 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 data sets (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 a test of the numerical methods (simplex and simulated annealing). Because these methods are indeed classical (Nelder & Mead 1965), 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 a test of systematic differences between the models rather than a test 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 that does not produce all the observables (astrometry, RVs, TTVs, _ecl, $| V{| }^{2}$, T₃, m_V, or I_λ) is impractical too, because we need as many 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 a such comparison in the future.

Finally, real observational data of real systems have their 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 subsequent system.

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 that are compatible within  respective approximations (e.g., on short timescales when the Keplerian model can be regarded as a useful approximation) with the following published works: the photometric model of Phoebe 1.0 (Prša & Zwitter 2006); the astrometric and speckle-interferometry model of Zasche & Wolf (2007); the RV disentangling by Korel (Hadrava 1997); the LitPro model for visibilities $| V{| }^{2}$ (Tallon-Bosc et al. 2008); the spectro-interferometric $| V{| }^{2}$ and T₃ (Nemravová et al. 2016); and the synthetic spectra fitting by Pyterpol (dtto).

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 χ² values) for Δt → 0.

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

Also, there is yet another discretization related to the WD code, or the surfaces of the eclipsing binary. For low numbers N_wd, one can see numerical artifacts on the light curve, as rectangular surface facets appear from behind the limb, or disappear. Again, it is worthwhile to check larger N_wd.

Quite often, we can expect one or more (m) mirror solutions (and 2^m combinations of them). A typical situation is that 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 Ω and Ω' = 180° − Ω are both admissible). Consequently, one may save some time when surveying the parameter space.

However, with the N-body model at hand it is worthwhile to check not only the total χ² but also individual contributions to χ² for all the mirror models. In particular,  ${\chi }_{\mathrm{ttv}}^{2}$ 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 data sets. If this is the unfortunate case, one may try to use weights w of individual χ²'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).

Radial-velocity measurements might be affected by zero-point offsets, which then lead to different systemic velocities γ 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 is atmospheric lines for which the relative RVs can be computed easily. If this is impossible, one should use the N-body model with great caution, because simply increasing σ_rv ≃ Δγ to get ${\chi }_{\mathrm{rv}}^{2}\simeq {N}_{\mathrm{rv}}$ is the wrong approach. 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 too much freedom, but if the dispersion relations can be considered stable from night to night, some calibration factors f_{rv k} of RVs—assigned to individual observatories or data sets—might actually be a better solution. In any case, such factors have to always be treated as additional free parameters of the N-body model.

Sometimes, RVs are derived in the Fourier domain by means of disentangling (e.g., by Korel; Hadrava 1995), with the advantage that one obtains 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 is a problem because we vary initial osculating orbital elements in the N-body model and they most likely will contradict the previous elements.

Note that 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 will obtain artificially small uncertainties σ_rv (and extremely large ${\chi }_{\mathrm{rv}}^{2}$). 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.

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_{rot j}, T_{eff j}, gravity $\mathrm{log}{g}_{j}$, and metallicity Z_j of the stellar components) are the same for all spectra, with the exception of RVs that 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 of its own when RVs are small (at conjunctions) and v_rot is large, which causes the lines to be totally blended. As a provisional solution, one may try to discard the lowest RVs that cause the problems, or not use RVs at all and instead fit synthetic spectra directly with the N-body model (${\chi }_{\mathrm{syn}}^{2}$), which is definitely a better approach, because RVs will be correctly tied to each other (see Figure 6).

Zoom In Zoom Out Reset image size

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

Contrary to closure phase arg T₃ 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 these suspicious observational data may be the only way to prevent the systematics to unrealistically shift the model. Using a low weight w_vis = 0.1 is not a satisfactory option. To this point, we always retain a data set identification for each single measurement that enables us to quickly perform a bootstrap testing.

A removal of quasi-periodic light oscillations that 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 around the minimum.

The observed light curve should thus be 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 }_{\mathrm{lc}}^{2}$; but there is no useful information, as we have no physical model for these oscillations (yet).

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 σ, because the former orbital elements are fixed, while the latter are only osculating initial conditions at t = T₀. Generally, all elements are time-dependent quantities, a₁(t), e₁(t), i₁(t), etc., whereas their oscillations are often larger than the 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.

It is also possible to run the N-body integrator separately, regardless of an observational time span, and study the long-term evolution and stability of stellar systems. We may wish to prefer those orbital solutions that 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 Δt_out is larger than half of the shortest orbital period, P₁/2.

In a modified version of the BS integrator (swift_bs_fp), we can use an online 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 (Šidlichovský & Nesvorný 1996) to extract proper elements. For N mutually interacting bodies, one can expect $2N$ eigenfrequencies of the system, which are usually denoted g_j and s_j. The corresponding amplitudes e_pj, $\sin \tfrac{1}{2}{i}_{{\rm{p}}j}$ can be considered approximate integrals of motion that only evolve on timescales longer than secular.

To conclude pessimistically, the above list of possible problems and systematics cannot be treated as complete, unfortunately.

By default, we assume a hierarchy of ((1 + 2) + 3) + 4, for which Jacobian orbital elements seem to be a suitable description. For a substantially different hierarchy, say, two pairs  (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_h2, y_h2, etc.), thus decreasing the dimensionality of the parameter space.

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 center of mass 1 + 2 is the reference point for the coordinates and velocities of the third body; the 1 + 2 + 3 center of mass is a suitable reference for the fourth 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° in the position angle—not due to an ambiguity, but simply because the reference body is different in our case. Similarly, a value of Ω from the literature may actually differ by 180°.