EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 666 (October 2022) A&A, 666 (2022) A24 Full HTML
Open Access
Issue
A&A
Volume 666, October 2022
Article Number A24
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202243596
Published online 30 September 2022

© M. Brož et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

HD 93206 is a complex system composed of nine components (usually denoted Aa1, Aa2, Ac1, Ac2, Ab, Ad, B, C, D; Mason et al. 2001). In the accompanying paper (Mayer et al. 2022) we used several observation-specific models to obtain the fundamental parameters of its subsystems, in particular the eclipsing binary Ac1+Ac2, the spectroscopic binary Aa1+Aa2, and their mutual long-period orbit (Aa1+Aa2)+(Ac1+Ac2). Together, the quadruple subsystem is called the QZ Car variable star. All these models were kinematical; they were based on the two-body problem with the light-time effect and the third light, to account at least for major perturbations.

In this work we focus on the same quadruple subsystem QZ Car, but we use a dynamical N-body model called Xitau1 to account for additional perturbations (in particular, all mutual, parametrized, post-Newtonian). Moreover, we account for all observational datasets at the same time, including astrometry, photometry, spectroscopy, and interferometry. This approach allows us to derive robust estimates of fundamental parameters. Our model was already used for multiple stellar systems (e.g., ξ Tau, Nemravová et al. 2016; Brož 2017) and asteroidal moon systems with highly irregular central bodies (e.g., (216), Marchis et al. 2021; Brož et al. 2021a, 2022).

2. Observational datasets

The radial velocity measurements, relative spectroscopy, light curves and corresponding eclipse-timing variations (ETVs), and astrometry are already described in Mayer et al. (2022). Here we focus on additional datasets, namely the spectral energy distribution (SED) and interferometry.

2.1. Spectral energy distribution

We used data from Vizier, which includes the standard Johnson system photometry (Ducati 2002), verified by the La Silla and Sutherland photometry (see Mayer et al. 2022), Hipparcos (Anderson & Francis 2012), Gaia3 (Gaia Collaboration 2020), 2MASS (Cutri et al. 2003), WISE (Cutri et al. 2012), MSX (Egan et al. 2003), and Akari (Ishihara et al. 2010). Altogether, the whole spectral range is 0.35–22 μm.

It was necessary to carry out the dereddening of these data. Unfortunately, the classic extinction maps of Green et al. (2019) give unrealistic absorption, and a special treatment is needed for QZ Car, which is too close to the galactic equator. To this end, we used the extinction maps of Lallement et al. (2019). For the assumed distance 2870 pc (Shull & Danforth 2019), the reddening extrapolated from 2100 pc is E(B − V) = (0.470 ± 0.060) mag, and corresponding absorption AV = 3.1 E(B − V) = (1.46 ± 0.19) mag (Fig. 1). Using the standard wavelength dependence of Schlafly & Finkbeiner (2011), we obtained the dereddened fluxes shown in Fig. 2.

thumbnail Fig. 1.

Reddening E(B − V) vs. distance for the direction towards QZ Car, according to Lallement et al. (2019). The mean (orange), dispersion (gray), and extrapolation (dashed) to the nominal distance 2870 pc are shown.

thumbnail Fig. 2.

Observed SED (blue) and dereddened SED (orange) for the distance 2870 pc. The dereddened fluxes Fλ approximately correspond to black-body radiation with the temperature 30 000 K (cyan). In the far-infrared, where absorption is negligible, a non-negligible (half order) excess is present.

Not surprisingly, the corrected fluxes exhibit a power-law slope, which corresponds to the Planck function (Bλ ∝ λ−4 for λ > 0.35 μm), for the effective temperature T ≃ 30 000 K. This is in agreement with the early-type spectral classification of QZ Car. Using a significantly smaller distance (e.g., 2000 pc) would result in a SED that does not correspond to early-type stars. We show in Sect. 3 that this distance and dereddening are indeed reliable. On the other hand, in the far-infrared (> 2.2 μm), where absorption is already negligible, a non-negligible excess is present. The monochromatic fluxes are a factor of 5 too large, possibly indicating cold circumstellar matter or blending, neither of which is accounted for in our modelling. Consequently, we should preferably fit the near-ultraviolet to near-infrared (NUV–NIR) spectral range.

2.2. Interferometry

In addition to the Sanchez-Bermudez et al. (2017) astrometric data, for the epochs June 17, 2016, and June 18, 2016, we used interferometric observations with the VLTI/GRAVITY instrument found in the ESO archive for the epochs March 14, 2017 and April 27, 2017. This is an important constraint for the (Ac1+Ac2)+(Aa1+Aa2) orbit. Apart from astrometry, we can directly include squared visibilities V2, closure phases argT3, and triple product amplitudes |T3| in Xitau.

We used the Esoreflex VLTI/GRAVITY pipeline (Freudling et al. 2013; GRAVITY Collaboration 2017) to reduce the data. We performed the standard visibility calibration, which is substantial (0.2) for baselines ≥4.5 × 107 cycles per baseline. We used the so-called ‘vfactor’ for the correction of V2. We did not use data with V2 + σV2 > 1.0. The new (u, v) coverage is shown in Fig. 3. The wide orbit ((Ac1+Ac2)+(Aa1+Aa2)), having the angular separation more than 30 mas, is clearly resolved. On the other hand, angular diameters of individual components (0.015–0.08 mas) cannot be resolved. As we discuss in Sect. 3, there were some problems with the |T3| calibration. Consequently, we also tried to reduce the data without ‘vfactor’.

thumbnail Fig. 3.

Coverage (u, v)≡B/λ (in cycles per baseline) of new interferometric observations (see Table 1). The corresponding squared visibility V2 is plotted in colour (see colour bar at right). The wide orbit ((Ac1+Ac2)+(Aa1+Aa2)), having the angular separation more than 30 mas, is clearly resolved.

As a check, we performed astrometry with the Litpro software (Tallon-Bosc et al. 2008). A simplified model of a binary (Ac+Aa) was assumed. We fitted (u, v) coordinates of Aa and the flux ratio FAa/FAc; other parameters were kept fixed. Even though there are some systematics, the fit converges without problems and the astrometric uncertainties are comparable to previous measurements (σu, v ≃ 0.010 mas; Table 1)2.

Table 1.

Astrometric positions of component (Ac1+Ac2) with respect to (Aa1+Aa2) from the literature and from new data.

We used no light curves in our modelling (only indirectly as ETVs). To our knowledge, there are no sufficiently high signal-to-noise ratio differential interferometric measurements (cf. Sanchez-Bermudez et al. 2017), which would otherwise constrain velocity-dependent visibilities across spectral lines.

3. Dynamical N-body model

In Xitau the dynamical model contains the following accelerations terms (Brož 2017; Brož et al. 2021a, 2022):

(1)

Here the terms are mutual gravitational interactions and the parametrized post-Newtonian (PPN); we do not activate oblateness, multipoles, or tides for QZ Car. We exchanged the PPN approximation for that of Standish & Williams (2006), which is more accurate. Hereafter, the notation conforms to the actual implementation:

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

Here ri denotes the position vector of body i, rij ≡ ri − rj, vi ≡ |ṙi|, c is the speed of light in a vacuum, and β = γ = 1 are the nominal PPN factors.

We perform a numerical integration of orbits with the Bulirsch-Stoer algorithm. It has an adaptive time step, with the relative precision set to 10−8. The output time step is set to 0.2 d, but all times of observations are also integrated and output precisely. The total number of free parameters is 38 (see Table 2). Observables are derived from coordinates, velocities, radiative quantities, and distance. Synthetic spectra, both relative and absolute, are computed during convergence with Pyterpol (Nemravová et al. 2016) from BSTAR and OSTAR grids (Lanz & Hubený 2007, 2003).

Table 2.

Best-fit parameters for the nominal (Sect. 3.2) and alternative (Sect. 3.3) models.

The observations and our model are compared by means of the χ2 metric

(14)

where the subscripts denote the respective datasets: radial velocity (RV), eclipse or transit timing variation (TTV), squared visibility (VIS), closure phase (CLO), triple product amplitude (T3), synthetic spectrum (SYN), spectral energy distribution (SED), and relative astrometry (SKY2). Recent improvements to our model also include metallicity fitting, a subplex algorithm (i.e. simplex on subspaces; Rowan 1990), and precise computation of the Roche potential from the volume-equivalent radius for optional light curve computations (using R(Ω) integrated as in Leahy & Leahy 2015 and inverted to Ω(R)).

As a preparatory task, we computed models with parameters as determined by the observation-specific models (Mayer et al. 2022). However, there were some caveats. First, orbital elements in our N-body model (including periods P1, P2, P3) are only osculating. When a dynamical model is different, these elements must be converged again. Whether or not fppn is included in the model also affects χ2; in other words, to obtain accurate values of the parameters it should be included. Moreover, we had to flip the orientation of the system several times to match some of the datasets (TTV, SKY2). Only then did we perform a convergence with the simplex algorithm (Nelder & Mead 1965; e.g., Fig. 4). This was performed multiple times using various initial conditions to avoid a false convergence and to solve some remaining systematics (RV, SYN).

thumbnail Fig. 4.

Example of convergence for the alternative model (see Table 2). All datasets (RV, TTV, VIS, CLO, T3, SYN, SED, SKY2) were taken into account, with alternative weights wrv = 100, wttv = 1000, wvis = 0.1, wclo = 0.1, wsyn = 0.1, wsky2 = 1000.0. Two of the stellar components (Ac2, Aa2) were forced to be on the main sequence, which results in tension (an increase in several χ2 contributions after 800 iteration). The nominal model does not exhibit such tension.

3.1. Survey of parameters

In order to understand how individual datasets constrain the model, and to find a global minimum of χ2, we perfomed a survey of the parameters. We converged 81 different models. Every convergence was initialized with a different combination of masses m1, m2 (Ac1, Ac2) in the range 15 to 50 MS, while m3 (Aa1) was set to msum − m1 − m2 − m4; nevertheless, all parameters were free during convergence. The maximum number of iterations was set to 1000. To speed up this computation, we used a set of predefined synthetic spectra. We assumed unit weights, with the exception of wsyn = 0.1. The radiative parameters of the fourth spectrally undetected component (Aa2) were constrained by Harmanec (1988) parametric relations m4(T4), log g4(T4) (i.e. assuming it is a normal main-sequence object). The component is too faint to significantly contribute to the total flux; it is not, however, negligible from the dynamical point of view.

The results are shown in Fig. 5. One can clearly see the forbidden regions where no solution can be found (χ2 remains high). For high m1, m2, this is especially due to the CLO and VIS datasets, which strongly constrain the relative luminosities of the (Ac1+Ac2) and (Aa1+Aa2) components. Moreover, there are strong correlations between parameters: negative for TTV, positive for RV. This is very useful because the best-fit model is just at the intersection. Finally, the χ2 contribution due to the SED dataset seems to be too flat. The reason is that the distance around 2800 pc corresponds to the SED, and all the models were converged. It does not mean that the SED is unimportant.

thumbnail Fig. 5.

Contributions to χ2 for a set of 81 best-fit models. Individual contributions (datasets) are shown in the panels (from top left): TTV, RV, SKY2, SED, SYN, VIS, CLO, T3. Every convergence was initialized with a different combination of masses m1, m2 (Ac1, Ac2) in the range 15–50 MS, while m3 (Aa1) was set to msum − m1 − m2 − m4. All parameters were free during convergence. The axes correspond to the masses m1 and m2; the colours correspond to χ2 (see also tiny numbers), with adapted colour scales: cyan best fits, blue good fits (< 1.2minχ2), orange poor fits (≥1.2minχ2). The factor was 3.0 for TTV, RV, SKY2. The forbidden regions can be clearly seen (e.g., high m1, m2 especially due to CLO), as well as correlations between the parameters (TTV −, RV +). The weighted very best fit for all the datasets is denoted by the red circle.

3.2. Nominal model

The best-fit parameters of our nominal model are listed in Tables 2 and 3. The observed and synthetic data are compared in Figs. 6–14.

thumbnail Fig. 6.

Eclipse timing variations (ETVs) for the nominal model (top; see Table 2). All synthetic minima (primary and secondary) are plotted on the x-axis (black), together with the observed minima (green) and their differences on the y-axis (red). In the N-body model these ETVs are small because they are suppressed by the actual motion of the eclipsing binary (Ac1+Ac2) about the centre of mass of all four components. For comparison, a simplified two-body model from Mayer et al. (2022), exhibiting large ETVs, is plotted (bottom; gray), with the light-time effect extracted from our N-body model (olive).

thumbnail Fig. 7.

Radial velocities (RVs) of components Ac1 (orange) and Ac2 (green). The observed values derived from the spectra (Mayer et al. 2022; error bars), synthetic values directly derived from our N-body model (dots), and residuals (red) are plotted. Large residuals are present close to the eclipses (phases 0.0, 0.5) because the individual lines (RVs) cannot be reliably separated, also for the component Ac2, which has very broad lines (hence uncertain RVs).

thumbnail Fig. 8.

Same as Fig. 7 for components Aa1 (blue) and Aa2 (magenta). The spectroscopic binary orbit (Aa1+Aa2) is highly eccentric. The broad range of the synthetic and the observed RVs is due to the motion on the wide orbit ((Ac1+Ac2)+(Aa1+Aa2)).

thumbnail Fig. 9.

Astrometry derived from previous interferometric measurements (see Table 1). Positions of the Aa1 components with respect to Ac1 are plotted (blue), as are Ac2 with respect to Ac1 (green; in the very centre). Tiny oscillations visible as thickness of lines are due to photocentre motions. We used positive signs of (u, v) for this plot. The orbit is constrained not only by this astrometry, but also by ETVs and RVs (see Figs. 6 and 8).

thumbnail Fig. 10.

Comparison of the rectified observed FEROS spectra (blue) with the synthetic spectra computed by our N-body model (orange). Residuals are also plotted (red). They are substantial for the Hδ line (see circles), where the continuum level is uncertain (between 405.1 and 412.6 nm). Nevertheless, most line features (including splits, blends, weak lines close to the continuum) are present in both observed and synthetic spectra, resulting in reliable RVs of individual stellar components (especially Ac1, Aa1; cf. Fig. 15).

thumbnail Fig. 11.

Comparison of the dereddened (blue) and the synthetic SED (red). The NUV–NIR region is reasonably described by our N-body model, although there is a substantial excess of the observed flux in the far-infrared region. This is possibly related to non-thermal processes in stellar atmospheres or to the presence of circumstellar matter, which is not included in our N-body model.

thumbnail Fig. 12.

Squared visibilities V2 for VLTI/GRAVITY observations (March 14, 2017 and April 27, 2017; blue), synthetic V2 (orange), and residuals (red). The sinusoidal dependence on the baseline B/λ essentially corresponds to a wide binary ((Ac1+Ac2)+(Aa1+Aa2)). Disks of individual stellar components cannot be resolved (at 0.08 mas, which correspond to a drop in V2 of only 0.01 at the longest baseline).

thumbnail Fig. 13.

Same as Fig. 12, but for the closure phase argT3. Again, the dependence on B/λ is similar to that of a binary.

thumbnail Fig. 14.

Same as Fig. 12, but for the triple product amplitude |T3|. The observed values are often close to or even exceed 1, which indicates a possible calibration problem at long baselines (B/λ > 4 × 107 cycles per baseline); synthetic values are always < 1. The positions of maxima and minima of |T3|(B/λ) are fitted perfectly.

thumbnail Fig. 15.

Synthetic spectra of individual components (Ac1, Ac2, Aa1, Aa2) weighted by their luminosities to show their contributions to the total spectrum. For comparison, the observed spectrum is plotted (error bars). The region of the He I 5016 and 5047 lines is plotted. The Ac2 component is fast-rotating and has broad lines; the Aa2 component is not luminous and its contribution is negligible.

Table 3.

Derived parameters for the nominal and alternative models.

First, the eclipse-timing variations (Fig. 6) exhibit much smaller amplitudes compared to those from Mayer et al. (2022) because they are suppressed by the actual motion of the eclipsing binary (Ac1+Ac2). The respective contribution is still larger than the number of data points nttv because the other datasets (RV) have to be fit at the same time; the result is a compromise.

Second, the RVs of the eclipsing binary Ac1+Ac2 (Fig. 7) show some systematics related to the Ac1 component at the phases 0.0 and 0.5 (i.e. primary and secondary eclipses) when the lines are blended and the RVs should be close to zero. Additional systematics are present elsewhere, for example between 0.6 and 0.7, which cannot be avoided due to other RV measurements with relatively low RV values. The reason might be hidden in a complex reduction and rectification procedure; one possible solution would be to prefer fitting the synthetic spectra, instead of deriving the RVs (see also Sect. 3.4.) For the Ac2 component, which is relatively faint and fast-rotating, a number of measurements are offset. Because these substantially contribute to , we decided to remove the RVs of Ac2 in the alternative model (see below).

Third, the RVs of the Aa1 component (Fig. 8) seem to be reliable, because it is relatively bright and slow-rotating. Even so, the Aa1 and Ac1 RVs are both clearly affected by the wide orbit ((Ac1+Ac2)+(Aa1+Aa2)). This was detected in Mayer et al. (2022) as variable γ velocities. The fourth component, Aa2, is too faint, but its predicted maximum RV should be of the order of ±300 km s−1.

Astrometry of the wide orbit (Ac+Aa binary; Fig. 9) provides a fundamental angular scale, which together with ETVs and RVs constrains the fundamental parameters. Even though there is some tension between periods from RVs and ETVs (Mayer et al. 2022), we regard the resulting period P3 = (14722 ± 50) d as a reasonable compromise.

The blue part (399–544 nm) of the eight rectified FEROS spectra (with high signal-to-noise ratio) were also fitted (Fig. 10). The contribution is still substantial, mostly because the level of continuum was not always correctly determined. The Hβ and Hγ lines were removed from the fitting procedure, due to strong emission features; it was impossible to be fit them by our model, which does not include any optically thin circumstellar matter (CSM). We can only fit weak emission lines (e.g., Fe III) present in the synthetic spectra of stellar atmospheres. On the other hand, Hδ does have features (asymmetries), which are reasonably fitted, although the depths of synthetic hydrogen lines are not matched precisely, which is again related to the level of the surrounding continuum; for example, between 405.1 to 412.6 nm there is no continuum, as verified by synthetic spectra. Nevertheless, the depths and the features of the He I, C II, O II, Mg II lines are fitted much better. Overall, the determination of the respective RVs is very precise (especially for the Ac1 and Ac2 components) because we fit all the lines at the same time in a process also known as cross-correlation. Our model also naturally couples all the spectra (by the orbits), which avoids line blends.

For the SED (Fig. 11), the fit is acceptable within ±0.2 mag in the NUV–NIR range, which is likely related to the uncertainties of dereddening (see Sect. 2.1). Moreover, 0.15 mag is the depth of the minima. The far-IR range was not fitted; the obseveved flux exhibits an excess, indicating a presence of the CSM. Nevertheless, the SED (together with the mass, RV, angular scale) determines the distance around 2800 pc.

A comparison of interferometric observations and our model is straightforward for the squared visibilities V2 (Fig. 12). The sinusoidal dependence on the baseline B/λ essentially corresponds to a wide binary ((Ac1+Ac2)+(Aa1+Aa2)). The scatter of the observed V2 is related to the Poisson noise; our synthetic V2 are smooth. Because V2 values can be very close to 1 (for a binary), it is important to include the observed values that are slightly higher or lower than 1 because only then will the mean value be close to 1. The corresponding precision of astrometry, when all VLTI/GRAVITY measurements are compressed to (u, v) coordinates, is of the order of 0.010 mas = 10 μas.

The same is true for the closure phases argT3 (Fig. 13). As in Sanchez-Bermudez et al. (2017), the dependence on B/λ is similar to that of a binary.

On contrary, the triple product amplitude (Fig. 14) exhibits some systematics, especially at long baselines (B/λ > 4 × 107 cycles per baseline); the observed values are often close to or even exceed 1, which seems to be a calibration problem. Nevertheless, the positions of the maxima and minima of |T3| as a function of B/λ are fitted perfectly, hence we retained the dataset in the nominal model (and we removed it in the alternative model; see Sect. 3.3).

3.3. Alternative model

Alternative models are useful for determining the true uncertainties of parameters. With this goal in mind, we modified a number of things, in particular the RVs of the Ac2 component were not used because they did not seem reliable. Instead, the Ac2 component was again constrained according to Harmanec (1988). In the case of the Ac1 RVs, we removed the outliers identified in the nominal model.

Similarly, |T3| was not used because it exhibited systematics at long baselines, while V2 at the same baselines was fitted precisely. Outliers were also removed from the VIS, CLO, and T3 datasets, as they often occur at the beginning of the burst. We thus assumed they were caused by a temporarily incorrect set-up.

The weights were set to extreme values, to enforce the fitting of datasets with a small number of points: wrv = 100, wttv = 1000, wvis = 0.1, wclo = 0.1, wsyn = 0.1, wsky2 = 1000.0. On the contrary, we relaxed the SED constraint by using wsed = 1.

The results are summarized in Table 2 and the figures can be found in Appendix A. The major differences between the models can be summarized as follows. Using Harmanec (1988) means forcing the Ac2 component to be on the main sequence. Since we needed a substantial mass ratio to explain the RVs of Ac1, the temperature of Ac2 increases (and its radius decreases). Putting less weight on the SED means that the distance is allowed to decrease (to 2500 pc) at the expense of systematic differences in the visible to NIR region. The ETVs are more centred with respect to the zero O − C value, including a group of measurements at JD ≃ 2458570. The RVs and astrometry of Aa1 indicate that the eccentricity of the wide orbit is still uncertain (log e3 = −0.28 vs. −0.54). A comparison of observed and synthetic spectra shows somewhat larger systematics for the He I lines. For interferometric quantities (V2, argT3), we detect only minor shifts in B/λ and offsets in amplitudes of the closure phase.

Even though the alternative model exhibits larger total χ2, it is a viable alternative and we can use it to estimate the uncertainties (see Table 2; last column).

3.4. A note on oblateness

The oblateness of components may also contribute to the precession of orbits. It can be characterized by the Love number k2, the rotation period P, and the body radius R (as in Fabrycky 2010), or alternatively by a series of multipoles, Cℓ, m, Sℓ, m, if the shape is more complex (as in Brož et al. 2021a). For binary stars, k2 is primarily determined by the density profiles and the Roche potential (Claret 2004). For soft bodies, with the polytropic index n = 3–4, the expected value is of the order of 10−2–10−3 (Yip & Leung 2017); it is even less for evolved stars with extended envelopes3.

According to our computations, the oblateness with k2 = 10−3 increases to 5660, but it can be compensated for by slightly adjusting the period P1. Only if k2 reaches ∼ 0.003 for the Ac1 component, which is probably the upper limit because log g1 ≲ 3.5, would the precession rate increase, hence as well as . Interestingly, the difference appears to be of the same order as the systematics of RVs mentioned in Sect. 3.2. Nevertheless, the precession is not as substantial; the observed RV curves do not exhibit a clear temporal trend, but rather a scatter (Mayer et al. 2022), and thus all the parameters (except P1) remain the same, within the uncertainties.

4. Discussion

Compared to previous works on QZ Car (Walker et al. 2017; Blackford et al. 2020), our new model indicates a higher total mass (137 vs. 112 MS) and a different distribution of masses within the binaries Aa1+Aa2 and Ac1+Ac2, preferring the mass ratios Aa1/Aa2 ≃8 and Ac1/Ac2 close to 1, respectively.

Given the photo-spectro-interferometric distance of (2800 ± 100) pc for the nominal model, QZ Car is likely related to the Carina Nebula (NGC 3372). However, when cluster distances were determined from the Gaia early Data Release 3 parallaxes (Shull et al. 2021; Göppl & Preibisch 2022), they turned out to be systematically smaller, possibly due to anomalous extinction, with RV ≡ AV/E(B − V) = 3.2–4.0. For Collinder 228, the median distance is 2470 pc and the angular size about 14′, which corresponds to the tangential size of only 10 pc. The radial size is not well constrained by parallaxes. This revised distance is in agreement with our alternative model of QZ Car.

Massive multiple stars in this region may shed light on possible progenitors of the η Carinæ event (Weigelt et al. 2016), which was suggested to be caused by mass transfer in triple systems, eventually leading to an explosive merger (Smith et al. 2018). Notably, the QZ Car quadruple system includes the eclipsing mass-transferring binary Ac1→Ac2, which will evolve substantially. According to our model, the mass ratio Ac2/Ac1 is relatively close to 1, which means that we are observing this binary at almost its closest distance. If, for example, the final mass of Ac1 decreases to 5 MS, the angular momentum conservation implies the final separation (as well as the envelope extent of Ac1; Paxton et al. 2015) of the order of 1000 RS (i.e. substantially less than the pericentre of (Aa1+Aa2)+(Ac1+Ac2) binary, a3(1 − e3) = 6300 RS). This does not lead to an imminent instability4. It is still unclear, what has been (and shall be) the role of loosely bound companions, Ab, Ad, B, C, D, imaged around QZ Car (Rainot et al. 2020) because we still lack their proper motions.

5. Conclusions

Using the Mayer et al. (2022) observation-specific models as initial conditions for convergence, we constructed a robust N-body model of the QZ Car quadruple system, fitting all three orbits and the stellar–radiative parameters at the same time. Because we used all types of observations, which are in a sense orthogonal, the model is well-behaved and does not exhibit strong correlations (or drifts). Independent constraints for individual parameters were only used when needed (e.g., for the luminosity–radius of the Aa2 component, which is too faint to be directly observable in the flux).

Our preferred model is the nominal one (see Table 2); the best-fit masses are m1 = 26.1 MS, m2 = 32.3 MS, m3 = 70.3 MS, and m4 = 8.8 MS, with uncertainties of the order of 2 MS, and the distance d = (2800 ± 100) pc. The alternative model, with the Ac2 component forced to be on the main sequence, exhibits some systematic differences, especially for the SED, which can be attributed to anomalous extinction with RV ∼ 3.4. We used it to determine the realistic uncertainties of parameters.

There are still some poorly constrained parameters, in particular the mutual inclinations of the orbits, because the longitude of nodes Ω1 is not constrained by eclipses. For simplicity, we assumed a co-planarity, but observations with the VLTI/GRAVITY instrument covering both short periods (P1, P2) could be used to constrain it; the astrometric uncertainties related to the Aa1 component are of the same order as the photocentre motions related to the eclipsing binary Ac1+Ac2.

As a future work we suggest including the observed light curve and strong emission lines (especially Hα and Hβ) in the model. However, this is generally difficult, given the possible presence of both optically thick and optically thin circumstellar matter (CSM). An appropriate treatment of the radiation transfer in the CSM is needed for this purpose (e.g., Brož et al. 2021b). Moreover, a presence of oscillatory signals in the TESS light curve also requires substantial model improvements (e.g., Conroy et al. 2020).

2

In Xitau, the (Ac1+Ac2) eclipsing binary is located in the centre, hence we use negative signs of (u, v) for the computations.

3

For comparison, an incompressible body has k2 = 0.75 (exact), and the Earth 0.295 (Lainey 2016).

4

Interestingly, the Aa1 component seems to be extremely massive (up to 70 MS), especially compared to the Aa2 component (less than 10 MS). This may be also related to a (putative) past mass transfer Aa2  →  Aa1.

Acknowledgments

M.B., P.H., and M.W. were supported by the Czech Science Foundation grant 19-01995S. We thank an anonymous referee for comments. We thank Michael Shull for a discussion about anomalous extinction.

References

  1. Anderson, E., & Francis, C. 2012, Astron. Lett., 38, 331 [Google Scholar]
  2. Blackford, M., Walker, S., Budding, E., et al. 2020, J. Am. Assoc. Var. Star Obs., 48, 3 [NASA ADS] [Google Scholar]
  3. Brož, M. 2017, ApJS, 230, 19 [Google Scholar]
  4. Brož, M., Marchis, F., Jorda, L., et al. 2021a, A&A, 653, A56 [Google Scholar]
  5. Brož, M., Mourard, D., Budaj, J., et al. 2021b, A&A, 645, A51 [EDP Sciences] [Google Scholar]
  6. Brož, M., Ďurech, J., Carry, B., et al. 2022, A&A, 657, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Claret, A. 2004, A&A, 424, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Conroy, K. E., Kochoska, A., Hey, D., et al. 2020, ApJS, 250, 34 [Google Scholar]
  9. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  10. Cutri, R. M., Wright, E. L., Conrow, T., et al. 2012, VizieR Online Data Catalog: II/311 [Google Scholar]
  11. Ducati, J. R. 2002, VizieR Online Data Catalog, Catalogue of Stellar Photometry in Johnson’s 11-color system [Google Scholar]
  12. Egan, M. P., Price, S. D., Kraemer, K. E., et al. 2003, VizieR Online Data Catalog: V/114 [Google Scholar]
  13. Fabrycky, D. C. 2010, in Exoplanets, ed. S. Seager, 217 [Google Scholar]
  14. Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Gaia Collaboration 2020, VizieR Online Data Catalog: I/350 [Google Scholar]
  16. Göppl, C., & Preibisch, T. 2022, A&A, 660, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93 [NASA ADS] [CrossRef] [Google Scholar]
  19. Harmanec, P. 1988, Bull. Astron. Inst. Czech., 39, 329 [NASA ADS] [Google Scholar]
  20. Ishihara, D., Onaka, T., Kataza, H., et al. 2010, A&A, 514, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Lainey, V. 2016, Celest. Mech. Dyn. Astron., 126, 145 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Lanz, T., & Hubený, I. 2003, ApJS, 146, 417 [NASA ADS] [CrossRef] [Google Scholar]
  24. Lanz, T., & Hubený, I. 2007, ApJS, 169, 83 [CrossRef] [Google Scholar]
  25. Leahy, D. A., & Leahy, J. C. 2015, Comput. Astrophys. Cosmol., 2, 4 [NASA ADS] [CrossRef] [Google Scholar]
  26. Marchis, F., Jorda, L., Vernazza, P., et al. 2021, A&A, 653, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466 [Google Scholar]
  28. Mayer, P., Harmanec, P., Zasche, P., et al. 2022, A&A, 666, A23 (Paper I) [Google Scholar]
  29. Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [Google Scholar]
  30. Nemravová, J. A., Harmanec, P., Brož, M., et al. 2016, A&A, 594, A55 [Google Scholar]
  31. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  32. Rainot, A., Reggiani, M., Sana, H., et al. 2020, A&A, 640, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Rowan, N. 1990, PhD Thesis, Univ. Texas Austin, USA [Google Scholar]
  34. Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
  35. Sanchez-Bermudez, J., Alberdi, A., Barbá, R., et al. 2017, ApJ, 845, 57 [NASA ADS] [CrossRef] [Google Scholar]
  36. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [Google Scholar]
  37. Shull, J. M., & Danforth, C. W. 2019, ApJ, 882, 180 [NASA ADS] [CrossRef] [Google Scholar]
  38. Shull, J. M., Darling, J., & Danforth, C. W. 2021, ApJ, 914, 18 [NASA ADS] [CrossRef] [Google Scholar]
  39. Smith, N., Andrews, J. E., Rest, A., et al. 2018, MNRAS, 480, 1466 [Google Scholar]
  40. Standish, E. M., & Williams, J. G. 2006, in Orbital Ephemerides of the Sun, Moon, and Planets, ed. P. Seidelmann (University Science Books) [Google Scholar]
  41. Tallon-Bosc, I., Tallon, M., Thiébaut, E., et al. 2008, in Optical and Infrared Interferometry, eds. M. Schöller, W. C. Danchi, & F. Delplancke, SPIE Conf. Ser., 7013, 70131J [CrossRef] [Google Scholar]
  42. Walker, W. S. G., Blackford, M., Butland, R., & Budding, E. 2017, MNRAS, 470, 2007 [NASA ADS] [CrossRef] [Google Scholar]
  43. Weigelt, G., Hofmann, K. H., Schertl, D., et al. 2016, A&A, 594, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Yip, K. L. S., & Leung, P. T. 2017, MNRAS, 472, 4965 [CrossRef] [Google Scholar]

Appendix A: Figures for the alternative model

For comparison, we show the corresponding figures for the alternative model (Fig. A.1 to Fig. A.8). Generally, they are similar to those of the nominal model, but they demonstrate how weights and additional constraints (for components Ac2 and Aa2) can change the best-fit solution.

thumbnail Fig. A.1.

Same as Fig. 6, but for the alternative model.

thumbnail Fig. A.2.

Same as Fig. 7, but for the alternative model.

thumbnail Fig. A.3.

Same as Fig. 8, but for the alternative model.

thumbnail Fig. A.4.

Same as Fig. 9, but for the alternative model.

thumbnail Fig. A.5.

Same as Fig. 10, but for the alternative model.

thumbnail Fig. A.6.

Same as Fig. 11, but for the alternative model.

thumbnail Fig. A.7.

Same as Fig. 12, but for the alternative model.

thumbnail Fig. A.8.

Same as Fig. 13, but for the alternative model.

All Tables

Table 1.

Astrometric positions of component (Ac1+Ac2) with respect to (Aa1+Aa2) from the literature and from new data.

In the text
Table 2.

Best-fit parameters for the nominal (Sect. 3.2) and alternative (Sect. 3.3) models.

In the text
Table 3.

Derived parameters for the nominal and alternative models.

In the text

All Figures

thumbnail Fig. 1.

Reddening E(B − V) vs. distance for the direction towards QZ Car, according to Lallement et al. (2019). The mean (orange), dispersion (gray), and extrapolation (dashed) to the nominal distance 2870 pc are shown.
In the text
thumbnail Fig. 2.

Observed SED (blue) and dereddened SED (orange) for the distance 2870 pc. The dereddened fluxes Fλ approximately correspond to black-body radiation with the temperature 30 000 K (cyan). In the far-infrared, where absorption is negligible, a non-negligible (half order) excess is present.
In the text
thumbnail Fig. 3.

Coverage (u, v)≡B/λ (in cycles per baseline) of new interferometric observations (see Table 1). The corresponding squared visibility V2 is plotted in colour (see colour bar at right). The wide orbit ((Ac1+Ac2)+(Aa1+Aa2)), having the angular separation more than 30 mas, is clearly resolved.
In the text
thumbnail Fig. 4.

Example of convergence for the alternative model (see Table 2). All datasets (RV, TTV, VIS, CLO, T3, SYN, SED, SKY2) were taken into account, with alternative weights wrv = 100, wttv = 1000, wvis = 0.1, wclo = 0.1, wsyn = 0.1, wsky2 = 1000.0. Two of the stellar components (Ac2, Aa2) were forced to be on the main sequence, which results in tension (an increase in several χ2 contributions after 800 iteration). The nominal model does not exhibit such tension.
In the text
thumbnail Fig. 5.

Contributions to χ2 for a set of 81 best-fit models. Individual contributions (datasets) are shown in the panels (from top left): TTV, RV, SKY2, SED, SYN, VIS, CLO, T3. Every convergence was initialized with a different combination of masses m1, m2 (Ac1, Ac2) in the range 15–50 MS, while m3 (Aa1) was set to msum − m1 − m2 − m4. All parameters were free during convergence. The axes correspond to the masses m1 and m2; the colours correspond to χ2 (see also tiny numbers), with adapted colour scales: cyan best fits, blue good fits (< 1.2minχ2), orange poor fits (≥1.2minχ2). The factor was 3.0 for TTV, RV, SKY2. The forbidden regions can be clearly seen (e.g., high m1, m2 especially due to CLO), as well as correlations between the parameters (TTV −, RV +). The weighted very best fit for all the datasets is denoted by the red circle.
In the text
thumbnail Fig. 6.

Eclipse timing variations (ETVs) for the nominal model (top; see Table 2). All synthetic minima (primary and secondary) are plotted on the x-axis (black), together with the observed minima (green) and their differences on the y-axis (red). In the N-body model these ETVs are small because they are suppressed by the actual motion of the eclipsing binary (Ac1+Ac2) about the centre of mass of all four components. For comparison, a simplified two-body model from Mayer et al. (2022), exhibiting large ETVs, is plotted (bottom; gray), with the light-time effect extracted from our N-body model (olive).
In the text
thumbnail Fig. 7.

Radial velocities (RVs) of components Ac1 (orange) and Ac2 (green). The observed values derived from the spectra (Mayer et al. 2022; error bars), synthetic values directly derived from our N-body model (dots), and residuals (red) are plotted. Large residuals are present close to the eclipses (phases 0.0, 0.5) because the individual lines (RVs) cannot be reliably separated, also for the component Ac2, which has very broad lines (hence uncertain RVs).
In the text
thumbnail Fig. 8.

Same as Fig. 7 for components Aa1 (blue) and Aa2 (magenta). The spectroscopic binary orbit (Aa1+Aa2) is highly eccentric. The broad range of the synthetic and the observed RVs is due to the motion on the wide orbit ((Ac1+Ac2)+(Aa1+Aa2)).
In the text
thumbnail Fig. 9.

Astrometry derived from previous interferometric measurements (see Table 1). Positions of the Aa1 components with respect to Ac1 are plotted (blue), as are Ac2 with respect to Ac1 (green; in the very centre). Tiny oscillations visible as thickness of lines are due to photocentre motions. We used positive signs of (u, v) for this plot. The orbit is constrained not only by this astrometry, but also by ETVs and RVs (see Figs. 6 and 8).
In the text
thumbnail Fig. 10.

Comparison of the rectified observed FEROS spectra (blue) with the synthetic spectra computed by our N-body model (orange). Residuals are also plotted (red). They are substantial for the Hδ line (see circles), where the continuum level is uncertain (between 405.1 and 412.6 nm). Nevertheless, most line features (including splits, blends, weak lines close to the continuum) are present in both observed and synthetic spectra, resulting in reliable RVs of individual stellar components (especially Ac1, Aa1; cf. Fig. 15).
In the text
thumbnail Fig. 11.

Comparison of the dereddened (blue) and the synthetic SED (red). The NUV–NIR region is reasonably described by our N-body model, although there is a substantial excess of the observed flux in the far-infrared region. This is possibly related to non-thermal processes in stellar atmospheres or to the presence of circumstellar matter, which is not included in our N-body model.
In the text
thumbnail Fig. 12.

Squared visibilities V2 for VLTI/GRAVITY observations (March 14, 2017 and April 27, 2017; blue), synthetic V2 (orange), and residuals (red). The sinusoidal dependence on the baseline B/λ essentially corresponds to a wide binary ((Ac1+Ac2)+(Aa1+Aa2)). Disks of individual stellar components cannot be resolved (at 0.08 mas, which correspond to a drop in V2 of only 0.01 at the longest baseline).
In the text
thumbnail Fig. 13.

Same as Fig. 12, but for the closure phase argT3. Again, the dependence on B/λ is similar to that of a binary.
In the text
thumbnail Fig. 14.

Same as Fig. 12, but for the triple product amplitude |T3|. The observed values are often close to or even exceed 1, which indicates a possible calibration problem at long baselines (B/λ > 4 × 107 cycles per baseline); synthetic values are always < 1. The positions of maxima and minima of |T3|(B/λ) are fitted perfectly.
In the text
thumbnail Fig. 15.

Synthetic spectra of individual components (Ac1, Ac2, Aa1, Aa2) weighted by their luminosities to show their contributions to the total spectrum. For comparison, the observed spectrum is plotted (error bars). The region of the He I 5016 and 5047 lines is plotted. The Ac2 component is fast-rotating and has broad lines; the Aa2 component is not luminous and its contribution is negligible.
In the text
thumbnail Fig. A.1.

Same as Fig. 6, but for the alternative model.
In the text
thumbnail Fig. A.2.

Same as Fig. 7, but for the alternative model.
In the text
thumbnail Fig. A.3.

Same as Fig. 8, but for the alternative model.
In the text
thumbnail Fig. A.4.

Same as Fig. 9, but for the alternative model.
In the text
thumbnail Fig. A.5.

Same as Fig. 10, but for the alternative model.
In the text
thumbnail Fig. A.6.

Same as Fig. 11, but for the alternative model.
In the text
thumbnail Fig. A.7.

Same as Fig. 12, but for the alternative model.
In the text
thumbnail Fig. A.8.

Same as Fig. 13, but for the alternative model.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.