GRAVITY Spectro-interferometric Study of the Massive Multiple Stellar System HD 93206 A

One of the most important observational clues to understand the formation of massive stars is multiplicity (see e.g., Apai et al. 2007; Zinnecker & Yorke 2007; Chini et al. 2012; Peter et al. 2012). Some plausible mechanisms to form massive multiples include (a) disk fragmentation (Bonnell & Bastien 1992; Monin et al. 2007; Krumholz et al. 2009), (b) the accretion of wider low-mass systems (Maeder & Behrend 2002; Bonnell & Bate 2005), (c) formation through failed mergers and stellar collisions in early dynamical interactions (Zinnecker & Bate 2002), and (d) disk-assisted capture (Bally & Zinnecker 2005). To discern which of the aforementioned models should be preferred, an analysis of the orbital properties in multiple systems (e.g., period distribution, mass ratios, orbital eccentricities, coplanarity, etc.) is required. This is not an easy task, because companions can cover spatial scales from a few 1/10 astronomical units (au) to 1000s of au, and contrast ratios of components can span several orders of magnitudes.

Sana et al. (2014) conducted the first statistical survey of massive multiple systems combining observations with high angular resolution facilities (in particular optical interferometry) with archival spectroscopic data. With a main sample of 96 O-stars, this study indicates that ∼90% of the massive stars are at least binaries and that ∼30% of them belong to a higher-degree multiple system. For example, the Trapezium system in Orion (Weigelt et al. 1999; Schertl et al. 2003; Kraus et al. 2009) is one of these systems that reveal the existence of stellar companions at multiple spatial scales. In fact, all of its four main components, like θ1 Ori B (Close et al. 2013), turn out to be embedded "micro-clusters" or binaries. However, it is not yet clear whether (i) these systems are stable on long timescales and (ii) how dynamical interactions with the companions can affect their evolution. Allen (2015) suggested that quintuplet systems, like θ1 Ori B, are prone to dissolve in around ∼100 crossing times, with some of the components merging or ejected (i.e., becoming the progenitors of runaway stars; see also Fujii & Portegies Zwart 2011).

Therefore, resolving the components of high-degree multiple systems, both in a statistical way and targeting individual systems, and monitoring their orbital motions are mandatory to answer the aforementioned questions. For this purpose, several high angular resolution techniques are required. Among them, infrared long-baseline interferometry is one of the most crucial techniques to characterize compact systems where the resolution provided by classical imaging techniques (e.g., speckle imaging or Adaptive Optics imaging) is limited. This is particularly true, now that the foremost K-band beam combiner GRAVITY has been successfully deployed at the Very Large Telescope Interferometer (VLTI). The new capabilities of GRAVITY not only offer us the possibility to observe high-mass multiples with milliarcsecond resolution, but also at an intermediate spectral resolution (R ∼ 4000) and higher sensitivity than previous interferometric instruments. In this study, we present the results of a GRAVITY spectro-interferometric study of HD 93206 A, an X-ray emitter that is known to be one of the most compact massive quadruple systems in the Galaxy. The object has been frequently investigated as a template for the formation of massive multiples.

1.1. The Source: HD 93206 A

HD 93206 (≡QZ Car) is a complex multiple system (see Figure 1) that is the brightest member of the open cluster Collinder 228, located in the Carina Nebula at a distance of 2.3 ± 0.1 kpc (Walborn 1995; Smith 2002, 2006). The known resolved components are listed in Table 1. Before proceeding, some clarifications are in order:

• 1.  
  There is some confusion in the literature regarding component nomenclature. Here, we follow that of the Washington Visual Double Catalog (WDS; Mason et al. 2001), which is the one that deals with the parent−child relationships most accurately.
• 2.  
  In the current WDS classification, the components of Aa–Ac form a quadruple system of two spectroscopic binaries with Ac being an eclipsing binary. In contrast, the nomenclature of Sana et al. (2014) includes both spectroscopic binaries in component Aa and identifies each member using numbers 1 to 4.
• 3.  
  The component E identified in the classification of Sana et al. (2014) corresponds to component D in the WDS.
• 4.  
  Component C could be a spurious detection. It has not been observed since 1934, and it is not detected as a 2MASS source.
• 5.  
  Components Ab, B, and D (and C if it is real) are significantly dimmer (${{\rm{\Delta }}}_{m}\gt 3.0$ mag) than Aa and Ac, which dominate the light output of the system and are the only ones bright enough to contribute to the integrated spectrum.

Zoom In Zoom Out Reset image size

Table 1.  Observed Stellar Pairs in HD 93206

Pair	Year	P.A.	Sep.	${\rm{\Delta }}m$	Reference
		(deg)	('')	(mag)
Aa, Ab	2014	324	1.00	3.9	Sana et al. (2014)
Aa, Ac	2016	328	0.03	0.4	This work
A, B	2012	276	7.07	5.8	Sana et al. (2014)
A, C	1934	93	8.80	⋯	Dawson & Aguilar (1937)
A, E	2012	303	2.58	7.2	Sana et al. (2014)

Note. Component nomenclature follows the WDS (see the text). The rest of the information comes from the listed references.

Download table as:  ASCIITypeset image

This paper deals with the two bright central binaries Aa and Ac (see Figure 1). Speckle observations with the CHARA speckle camera (Mason et al. 1998) and direct imaging with the Fine Guidance Sensor of the Hubble Space Telescope (Nelan et al. 2004) could not resolve Aa–Ac, setting an upper limit of $\rho \sim 15.2$ mas (∼35 au) for their separation. However, H-band interferometric data taken with PIONIER (Precision Integrated-Optics Near-infrared Imaging ExpeRiment) at the VLTI and K-band sparse aperture masking (SAM) observations obtained with the near-infrared camera NACO at the VLT resolved the separation between the two binaries Aa and Ac, finding a projected angular distance of $\rho \sim 26$ mas (64.4 au in projection) and with brightness ratios of ${\rm{\Delta }}H=0.42$ mag and ${\rm{\Delta }}{K}_{{\rm{S}}}=1.18$ mag (Sana et al. 2014).

The spectroscopic analysis of the system (Leung et al. 1979; Morrison & Conti 1980; Stickland & Lloyd 2000; Mayer et al. 2001) indicates the existence of two periods of 20.7 and 6 days for components Aa and Ac, respectively. Therefore, from the spectroscopic point of view, HD 93206 Aa–Ac is a quadruple system that can be subdivided into two short-period binary systems, Aa1, Aa2, and Ac1, Ac2, in a long-period orbit of ∼50 years around each other. This configuration is different from other hierarchical quadruple systems (such as HD 17505; Hillwig et al. 2006; Sota et al. 2011) composed of four O-stars in which the components are organized in a close pair (with a period of days) orbiting a third star with a period of years and a fourth star orbiting around the other three with a period possibly measured in millennia. The intriguing architecture of the system highlights the importance of characterizing the orbital motion of the components in HD 93206 A with the final goal of confronting the different formation scenarios with its dynamical properties. For example, core collapse massive star formation models (e.g., Krumholz et al. 2009) assume that companions have coplanar orbital planes, while competitive accretion models (e.g., Bonnell & Bate 2006) predict non-coplanar orbits.

Spectroscopic analysis also suggests that Ac1, Ac2 is in fact a semi-detached eclipsing binary, and Aa1, Aa2 is identified as a long-period spectroscopic binary. The individual spectral classification of the stellar components in HD 93206 Aa–Ac supports Aa1 being an O-type supergiant with spectral type O9.7 Iab or O9.5 I (the classification for the combined spectral type is O9.7 Ibn; Sota et al. 2014). Aa2 is undetected in the spectrum but indirect evidence points toward it being several times less massive than Aa1 and with a spectral type in the vicinity of B2 V. Ac1 appears to be a giant or supergiant (O8 III or B0 Ib). Although the primary Ac1 is the brightest component, it seems to be less massive than the secondary Ac2 (possibly O9 V, but its signature on the integrated spectrum is weak), which suggests case-B mass transfer within this binary. The whole system has an estimated total mass of ∼90 M_⊙ (Leung et al. 1979; Morrison & Conti 1980; Mayer et al. 2001).

X-ray observations of the quadruple Aa–Ac (Broos et al. 2011; Townsley et al. 2011) revealed an excess in X-ray emission (average ∼7 × 10⁻¹³ erg s⁻¹ cm⁻²) over the expected cumulative bolometric luminosity. A three-component plasma model fitted to the observed X-ray spectra by Parkin et al. (2011) supports the presence of shock-heated plasma from wind–wind collision shocks. Those authors suggest that the system Ac is the main X-ray emitter. However, the total observed X-ray flux might be the sum of the emission of the individual components and the mutual wind–wind collision region between the two binaries. If this hypothesis is correct, we expect to observe infrared shock tracers (like Brγ) to be formed in an extended region where the colliding winds of the two binary systems interact.

The paper is outlined as follows. Section 2 presents our GRAVITY/VLTI observations and data reduction. In Section 3, the analysis of the interferometric observables and of the source spectrum are described. In Section 4, our results are discussed and, finally, in Section 5, the conclusions are presented.

2.1. GRAVITY Observations

HD 93206 A was observed for a total of two hours with GRAVITY (Eisenhauer et al. 2008, 2011; GRAVITY Collaboration et al. 2017) during the ESO Science Verification run of the instrument.⁷ Snapshot observations were obtained at four different hour angles, two of them during 2016 June 17 and two more during 2016 June 18. The observations were carried out using the highest spectral resolution (R ∼ 4000) of the interferometer in the K-band (1.990–2.450 μm), together with the combined polarization and single-field modes of the instrument. With this configuration, GRAVITY equally splits the incoming light of the science target between the fringe tracker and the science beam combiners to produce interference fringes simultaneously in each of them. While the science beam combiner disperses the light at the desired spectral resolution, the fringe tracker beam combiner always works with a low-spectral resolution of R ∼ 22 (Gillessen et al. 2010) but at a high-frequency sampling (∼1 kHz). This allows for the atmospheric piston to be corrected, stabilizing the fringes of the science beam combiner for up to several tens of seconds.

All of the HD 93206 A data were recorded in the A0-G1-J2-K0 configuration of the 1.8 m Auxiliary Telescopes (ATs) of the VLTI. This configuration provides a maximum projected baseline of ∼120 m (J2-A0) and a minimum one of ∼39 m (K0-J2) for the target coordinates. These baseline lengths correspond to a maximum angular resolution (θ = λ/2B) of $\theta \sim 1.89$ mas and a minimum one of $\theta \sim 5.8$ mas, at a central wavelength of ${\lambda }_{0}=2.2$ μm. Figure 2 displays the science target u − v coverage obtained with the aforementioned GRAVITY observations.

Zoom In Zoom Out Reset image size

The interferometric observables (squared visibilities, closure phases, and differential phases) as well as the source's spectrum were obtained using version 0.8.4 of the instrument's data reduction software⁸ provided by ESO (Lapeyrere et al. 2014); it uses a series of esorex⁹ recipes to estimate the raw observables and to calibrate them. The data reduction process converts the science camera frames (where the interference fringes are recorded using the ABCD sampling method; Colavita et al. 1999) into a set of 1D spectra from which the observables are estimated.

First, the science detector is characterized by measuring its bias, pixel gain, bad pixel and wavelength maps. The wavelength maps of the Fringe Tracker and Science beam combiners are calibrated using the GRAVITY Calibration Unit (Blind et al. 2014). For the Science Beam Combiner, the fiber wavelength scale is obtained through observations of the internal halogen lamp with the fiber differential delay lines in close loop and fringe tracking. At the same time, the calibration unit delay line positions are modulated, while their positions are monitored using the internal laser diode (${\lambda }_{0}=1.908$ μm). The measured fringe phase shifts are converted to wavelength knowing the introduced optical path delay (OPD) offsets, with a wavelength precision of ${\rm{\Delta }}\lambda =2$ nm. An argon spectral calibration lamp (which provides 10 lines in the GRAVITY wavelength range) is used to derive the vacuum wavelength scale and correct the wavelength map of the fiber dispersion, achieving an absolute wavelength calibration equivalent to one-half of the spectral element at the highest spectral resolution of the instrument (0.1 nm, equivalent to a relative uncertainty of $5\times {10}^{-5}$ in high spectral resolution mode)

The transfer function of the science beam combiner is calibrated by using the P2VM algorithm (Tatulli et al. 2007; Lacour et al. 2008). It first requires the V2PM matrix to be determined, as defined in Lapeyrere et al. (2014), for each spectral bin. The internal halogen lights are used to calibrate this matrix in three steps: first, for the power transmission of each of the telescopes; second, by computing the visibility amplitude and phase for each baseline (i.e., internal contrast and phase) produced by an induced OPD modulation; and third, by computing the closure phases. Finally, the inverse of the V2PM matrix is used to extract the interferometric observables from on-sky observations. The P2VM calibration files are generated from the raw data with the esorex recipe gravity_p2vm. Once they are created, the gravity_vis recipe can be used to extract the total flux per telescope and correlated flux per baseline (i.e., the raw observables). The recipe gravity_vis includes a bootstrapping method to compute the complex visibility and bispectrum errors, as well as a visibility-loss correction factor due to the integration time difference between the fringe tracker frames and the science ones.

To calibrate the interferometric observables, interleaved observations of the science target and point-like sources were performed. Table 2 displays the observational setup used for the science target together with the different calibrators. The interferometric calibration was performed using the gravity_viscal recipe. This routine computes the instrumental transfer function by correcting the observed calibrator visibilities with the theoretical ones according to the estimated angular size of the calibrator. The algorithm then groups calibrators with the same observational setup as the science target and interpolates their transfer functions to the acquisition time of the science target observations. Finally, the recipe corrects the target raw visibilities for the estimated transfer function, delivering the calibrated observables.

Table 2.  GRAVITY/VLTI Science Target and Calibrator Data Sets

Date	Target	Typea	${K}_{\mathrm{mag}}$	DITb	NDITc	No. PAd	Airmasse	Seeinge	UD sizef
2016 Jun 17	HIP 50644	CAL	4.25	30	10	2	1.49, 1.54	0.45, 0.51	0.853
	HD 93206 A	SCI	5.25	30	10	2	1.54, 1.67	0.57, 0.54	⋯
	HD 94776g	CAL	3.40	10	26	2	1.73, 1.79	0.80, 0.65	0.917
	HD 149835	CAL	4.90	30	10	2	1.17, 1.19	0.34, 0.64	0.49
	HD 188385	CAL	5.99	30	10	2	1.45, 1.51	0.47, 0.43	0.224
2016 Jun 18	HD 93206 A	SCI	5.25	30	10	2	1.40, 1.43	1.0, 0.77	⋯
	HD 94776g	CAL	3.40	10	26	2	1.45, 1.47	0.57, 0.53	0.917
	HD 166521	CAL	4.52	30	10	2	1.06, 1.07	0.48, 0.48	0.578
	HD 164031	CAL	4.11	30	10	2	1.31, 1.37	0.59, 0.88	0.707
	HD 8315	CAL	4.04	30	10	1	1.18	0.78	0.749

Notes.

^aThis parameter indicates the purpose of the target observed. CAL means that the source corresponds to a calibrator star and SCI means that the source corresponds to our science target. ^bDetector integration time in seconds. ^cNumber of frames to estimate the complex visibilities per data cube. ^dNumber of observed position angles per night. ^eAirmass and seeing per target data set. ^fK-band estimated uniform disk size in milliarcseconds obtained from SearchCal (Bonneau et al. 2006, 2011). ^gK0III star used to calibrate the HD 93206 A spectrum.

Download table as:  ASCIITypeset image

Before analyzing the data, those V² points with a signal-to-noise ratio (S/N) ≤ 5 and closure phases with ${\sigma }_{{cp}}\geqslant 40^\circ$ were flagged. Figures 11 and 12 in the Appendix display the distributions of the V² S/N and of the errors in the closure phases after flagging. From these distributions, it is observed that three baselines have a considerably small S/N, particularly in the second data set where only a few points remain after flagging. These baselines are connected with the K0 station, which shows a poor transmission throughput at the time of the Science Verification observations. The first data set has the highest quality mainly because of the excellent weather conditions under which it was taken (see Table 2).

2.2. FIRE Observations

2.3. FEROS Observations

3.1. Modeling the Interferometric Observables

3.1.1. MCMC Modeling

In the calibrated V² and closure phases, the target presents the prototypical cosine signature of a resolved binary system. Figure 3 displays one of the GRAVITY data sets (MJD:57557.046) where the binary signature is clearly visible. As expected, with the current GRAVITY observations, we could only map the separation between the two components Aa and Ac, but not resolve the individual components of these binaries (which have angular separations considerably less than 1 mas). A geometrical model fitting of a binary composed of two point-like sources was applied to the observables. The expression used to estimate the complex visibilities, $V(u,v)$, is the following:

$$\begin{eqnarray}&&V(u,v)=\displaystyle \frac{1+{F}_{\mathrm{Ac}}/{F}_{\mathrm{Aa}}G(\zeta ){e}^{-2\pi j({\rm{\Delta }}{\theta }_{{\rm{x}}}u+{\rm{\Delta }}{\theta }_{{\rm{y}}}v)}}{1+{F}_{\mathrm{Ac}}/{F}_{\mathrm{Aa}}},\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}&&G(\zeta )=\left|\displaystyle \frac{\sin \zeta }{\zeta }\right|\,\,\,\mathrm{with}\,\,\,\zeta =\displaystyle \frac{\pi ({\rm{\Delta }}{\theta }_{{\rm{x}}}u+{\rm{\Delta }}{\theta }_{{\rm{y}}}v)}{R},\end{eqnarray} \tag{ 2 }$$

where ${F}_{\mathrm{Ac}}/{F}_{\mathrm{Aa}}$ is the flux ratio between the components of the binary, (u, v) are the components of the spatial frequency sampled per each observed visibility point (u = B_x/λ and v = B_y/λ), and (${\rm{\Delta }}{\theta }_{{\rm{x}}}$, ${\rm{\Delta }}{\theta }_{{\rm{y}}}$) are the east–west and north–south components of the projected angular separation between the primary and secondary components of the binary. $G(\zeta )$ is a correction factor due to the bandwidth smearing of the finite bandpass used, and R = λ/Δλ is the spectral resolution of the interferometer (in this case R ∼ 4000). The total flux of the binary is normalized to unity such that F_Aa+F_Ac = 1.

Zoom In Zoom Out Reset image size

V² and closure phases across the entire bandpass were simultaneously used to estimate the average ${F}_{\mathrm{Ac}}/{F}_{\mathrm{Aa}}$ and (${\rm{\Delta }}{\theta }_{x}$, ${\rm{\Delta }}{\theta }_{y}$) of the binary. The fitting process was performed using a dedicated Markov chain Monte Carlo (MCMC) routine based on the Python software emcee (Foreman-mackey et al. 2013). First, an approximation of the best-fit parameters was obtained using the nonlinear least-squares minimization algorithm inside the Python library scipy.optimize.¹⁰ To account for the standard deviations and correlations of the model parameters, we explore the solution space around the best-fit values obtained with the least-squares method; 300 random points with a linear distribution between ±10σ were generated for each of the parameters. For every random point, 700 iterations were performed using the MCMC algorithm to maximize the posterior probability of the model. The ${\chi }^{2}$ of the best-fit model is ${\chi }^{2}=3.93$, while the parameter values and 3σ uncertainties of the best-fit model are reported in Table 3. Figure 4 shows the posterior probability distributions of the parameters, their correlations, as well as their individual marginalized distributions in 1D histograms along the diagonal of the plot.

Zoom In Zoom Out Reset image size

Table 3.  Best-fit Geometrical Model with the Different Algorithms Used in This Work

	MCMC Fit		LitPro Fit		CANDID Fit		Average Fit
Parameter	Value	±3σ	Value	±3σ	Value	±3σ	Value	±3σ
${F}_{\mathrm{Ac}}/{F}_{\mathrm{Aa}}$	0.710	±0.002	0.67	±0.11	0.832	±0.004	0.740	±0.14
${\rm{\Delta }}{\theta }_{x}$ [mas]	−15.965	±0.0013	−15.91	±0.020	−15.942	±0.0054	−15.940	±0.045
${\rm{\Delta }}{\theta }_{y}$ [mas]	25.593	±0.0034	25.645	±0.024	25.591	±0.005	25.609	±0.054

Download table as:  ASCIITypeset image

3.1.2. LitPro Modeling

In addition to our MCMC modeling, we also performed the analysis of the binary using LitPro¹¹ (Tallon-Bosc et al. 2008). This is a code that uses a gradient descent method to perform the optimization; it is based on a Levenberg–Marquardt algorithm combined with a Trust Region method. The software allowed us to use several pre-defined geometrical components to represent the interferometric data. In this case, we use a setup similar to the MCMC model (a binary composed of two point-like sources). There are two main differences between the two methods: (i) the first one is the lack of a correction in LitPro due to bandwidth smearing, and (ii) the second one is the flagging method applied by LitPro, which masks ${V}^{2}\lt 3{\sigma }_{{V}^{2}}$ or ${V}^{2}\gt 1+3{\sigma }_{{V}^{2}}$. Table 3 displays the best-fit model obtained with this software with an achieved ${\chi }^{2}=3.89$. Despite the differences in the algorithms and their setup, the parameters obtained with MCMC and LitPro are in very good agreement.

3.1.3. CANDID Modeling and Detection Limits

We also used CANDID¹² (Gallenne et al. 2015) to analyze the interferometric data of HD 93206 A. This Python code uses a Levenberg–Marquardt algorithm to fit for a binary system performing a systematic exploration over a given grid. The model of the binary represents a spatially resolved primary component (to which a uniform disk is fitted), as well as the bandwidth smearing of the visibilities (see Lachaume & Berger 2013). In this case, we considered the two components of the binary to be unresolved, and the size of the binary component was set to zero. For the grid search, we selected a step of p = 1.89 mas, which is equivalent to λ/2B_max. In addition to the binary fitting, CANDID also allows an estimation of the 3σ detection limit for the companions at different separations from the primary to be obtained. The code includes the algorithms described in Absil et al. (2011) and Gallenne et al. (2015) to look for high-contrast companions in the searched grid.

Table 3 displays the best-fit parameters obtained with CANDID. Figure 5 displays the map of the ${\chi }^{2}$ minima of the binary model with the position of the secondary marked with a red cross. In the case of our data, the secondary was detected with the highest confidence value allowed by CANDID (50σ). The minimum ${\chi }^{2}$ achieved by the searching algorithm is ${\chi }^{2}=2.46$. The flux ratio and position of the secondary found with CANDID agree with the values obtained with the other two methods previously described (see Table 3).

Zoom In Zoom Out Reset image size

CANDID was also used to estimate the maximum 3σ magnitude difference, ΔMag_3σ, to detect a companion for a separation range between 3 and 50 mas using both detection methods available in the code. Here, we only used one of the GRAVITY data sets (MJD:57557.046) to properly compare the results with the 2012 PIONIER data set reported by Sana et al. (2014). The published PIONIER astrometric position and flux ratio were adopted to subtract the secondary component from the data and search for faint companions in the grid with a step of p = 1.4 mas. Figure 6 displays the resulting detection limits. The blue and red solid lines indicate the average detection limits between the two detection methods available in CANDID. These results support the maximum detectable magnitude difference in the PIONIER data being ∼3 mag, while for GRAVITY it goes up to ∼5 mag.

Zoom In Zoom Out Reset image size

As evident in Table 3, the model parameters obtained with the three different methods agree quite well. However, the standard deviation obtained for each of the algorithms appears to be underestimated. This is mainly caused by slight deviations in the expressions used for the optimization, as well as from the flagging methods applied in each of the methods. To obtain a more robust estimation of the best-fit parameters and their uncertainties, we computed the mean and the standard deviation among the best-fit models determined with the different optimization methods. These estimates are also reported in columns 8 and 9 in Table 3, and those are the values adopted for the following analysis in this work. Figure 3 displays (for clarity) one of the GRAVITY data sets together with the best-fit average model overplotted. Notice how the trend in V² and CPs are well reproduced by the best-fit model.

3.2. Extracting the Source Spectrum

Together with the interferometric observables, GRAVITY delivers the spectrum of the source for each of the four telescopes that form the interferometer. The esorex data reduction software delivers the spectrum flattened by the instrumental transfer function. To correct for atmospheric effects, it is necessary to use a spectroscopic calibrator. Here, the K0III star HD 94776 was used (see Table 2). The spectral type of this star does not exhibit strong absorption or emission lines between 2.0 and 2.29 μm. However, it shows a CO band head in absorption from 2.29 μm onward. Therefore, wavelengths larger than 2.29 μm were discarded from the spectral analysis.

The calibrated spectrum of HD 93206 A was obtained by dividing each of the four raw-science spectra per data set (telescope) by its corresponding raw-calibrator spectrum. Every science data set was corrected by each of the calibrator observations per night, obtaining a total of 32 samples of the source spectrum. After the ratio was computed, the different calibrated spectra were normalized with a third degree polynomial fitted to the continuum. Due to the presence of outliers in the different spectra (associated with bad pixels) and to the difference in the S/N caused by different telescope transmissions, a Principal Component Analysis (PCA) was used to obtain a median spectrum of the source.

The PCA algorithm (see Jolliffe 1986) filters out the high frequencies in the data that correspond mostly to noise, allowing us to keep the components of the signal that maximize the variance (i.e., the direction with the maximum signal) between the different samples of the HD 93206 A spectrum. The PCA algorithm (Ivezic et al. 2014) consists of the following steps:

• 1.  
  We center each of the spectra in the set of data $\{{x}_{i}\}$ by subtracting the mean spectrum.
• 2.  
  A matrix ${\boldsymbol{X}}$ of N × K dimensions was built with the previously subtracted spectra $\{{x}_{i}^{\mathrm{sub}}\}$. Here, N corresponds to the number of samples of the HD 93206 A spectrum (i.e., N = 32) and K to the number of dimensions in the spectrum (i.e., to the spectral bins in the data).
• 3.  
  PCA aims to project ${\boldsymbol{X}}$ into space ${\boldsymbol{Y}}={\boldsymbol{XR}}$, where ${\boldsymbol{R}}$ is aligned to the direction of maximal variance. This projection is done through a single value decomposition (SVD) of ${\boldsymbol{X}}$ itself in its eigenvalues and eigenvectors.
• 4.  
  The reconstruction of each individual spectrum, ${x}_{i}$, is given by

  $$\begin{eqnarray}&&{{\boldsymbol{x}}}_{i}(k)={\boldsymbol{\mu }}(k)+\displaystyle \sum _{j}^{r\lt R}{\theta }_{{ij}}{{\boldsymbol{e}}}_{j}(k),\end{eqnarray} \tag{ 3 }$$

  where the indices k correspond to the different wavelength bins, i to the number of the input spectrum, and j to the number of eigenvectors used. ${\boldsymbol{\mu }}$ corresponds to the mean of the original spectra per wavelength bin. R is the total number of eigenvectors ${\boldsymbol{e}}(k)$. In this case, the summation was performed only up to the coefficient r which corresponds to a cumulative variance of 50% in the data. The coefficients, ${\theta }_{{ij}}$, are given by

  $$\begin{eqnarray}&&{\theta }_{{ij}}=\displaystyle \sum _{k}{{\boldsymbol{e}}}_{j}(k){{\boldsymbol{x}}}_{i}^{\mathrm{sub}}.\end{eqnarray} \tag{ 4 }$$

Figure 7 displays a normalized median of the HD 93206 A spectrum after the PCA analysis. The displayed error bars correspond to the standard deviation of the median. Although no prominent spectral features are observed, some of them, like Brγ and He i, are identified. For comparison, a normalized spectrum of the source taken with the FIRE spectrograph, at the 6 m Magellan telescope, is also plotted. Similar spectral features are distinguished in both spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the geometrical model, a separation between the components in the outer binary of ρ = 30.16 ± 0.02 mas and a position angle of θ = 32809 ± 003 were derived. These estimates agree with the values reported by Sana et al. (2014) using the PIONIER (ρ = 25.76 ± 0.54 mas; θ = 33131 ± 145) and SAM (ρ = 28.05 ± 5.41 mas; θ = 33102 ± 961) observations taken in 2012, but with a precision between one and three orders of magnitude better, respectively. The SAM data derived ${\rm{\Delta }}{K}_{s}=1.18\pm 0.20$ mag; however, our GRAVITY observations suggest ΔK = 0.33 ± 0.02 mag, a value that is closer to the ΔH = 0.42 ± 0.18 mag derived with the PIONIER observations. The angular separation between the components of the binary is at the limit of the SAM angular resolution (${\theta }_{K}\sim 35$ mas, for 8 m telescopes). Thus, this could explain the poor performance of such a technique to reliably recover the ${\rm{\Delta }}K$ and separation of the system Aa–Ac. In addition, it highlights the importance of long-baseline interferometric observations, like our GRAVITY data, to cover the relevant angular scales of systems like HD 93206 A.

With the position angles derived in 2012 and in this work, a projected northwest motion in the plane of the sky of component Ac relative to Aa is detected. Figure 8 displays a BSMEM (Buscher 1994; Baron & Young 2008) reconstructed image from our GRAVITY data together with the 2012 PIONIER and SAM positions marked on it. With only two astrometric points, it is not possible to derive a reliable orbital solution of the outer binary Aa–Ac (see Figure 8). Nevertheless, (i) assuming a total mass of $90\,{M}_{\odot }$ (Leung et al. 1979; Morrison & Conti 1980), (ii) a projected circular orbit with a mean separation at the aphelion of ∼30 mas, (iii) and a parallax to the Carina Nebula $\epsilon =0.43\pm 0.02$ (Smith 2002), a rough orbital period of P = 61 ± 4 a is inferred. However, this estimate highly depends on the orbital elements of the system (e.g., eccentricity, inclination, mass, etc.). For example, taking the previous constraints into account, a system projected in the plane of the sky with eccentric orbits e = 0.5 and e = 0.9 will result in P = 33 ± 2 a and P = 23 ± 1.6 a, respectively. These, still, poor constrains on the orbit support the necessity of carrying out a proper GRAVITY monitoring program over, at least, the next five years, which, together with the previous 2012 and 2016 data, would allow us to obtain more accurate estimates of the orbital motion of the system. The current Gaia TGAS (Gaia Collaboration et al. 2016a, 2016b) lists for HD 93206 a highly uncertain parallax of $\epsilon =1.76\pm 0.62$ mas, which includes, e.g., the distance to the Carina nebula within its 2σ uncertainty. Once more a precise Gaia parallax measurement becomes available, it (combined with our new astrometric data) will further improve the constraints on the mass of the stellar components.

Zoom In Zoom Out Reset image size

The general spectral signature of HD 93206 A obtained with GRAVITY is consistent with the one taken with FIRE. In both spectra, the Brγ (2.1661 μm) shock tracer and He i (2.058 and 2.112 μm) photospheric lines can be identified. The observed line profiles resemble those of late O- and early B-stars (see, e.g., Hanson et al. 1996; Ghez et al. 2003; Tanner et al. 2006). This result confirms the presence of O9.7 I and O8 III stars in HD 93206 A. However, some differences between the FIRE and GRAVITY spectra are observed. The He i line at 2.112 μm in the FIRE spectrum appears to be an absorption feature with a well-defined peak and a wing shifted toward the red part of the spectrum. In contrast, the GRAVITY spectrum shows an absorption "W-shaped" profile. One of the most plausible causes of these line changes is the difference of the orbital phase in the compact binaries between both spectra (see, e.g., Figures 1 and 2 in Taylor et al. 2011). Variability of the spectral features in the visible part of the spectrum of HD 93206 A was also reported by Morrison & Conti (1979) and Mayer et al. (2001), associating it with the double tight-binary nature of the source.

To test this scenario, Figure 9 shows normalized spectra in the Hα line (0.656 μm) of HD 93206 A obtained with FEROS (see Section 2.2). Considering the ephemeris published by Mayer et al. (2001) for Aa+Ac, we determine the orbital phases corresponding to each of the spectroscopic binaries. The different phases are labeled to the left and right of each spectrum, respectively. At first glance, the variations of the profile on a day-to-day basis are very noticeable. The 2009 series of spectra (the four spectra from the bottom to the top in Figure 9) were obtained after the periastron passage of component Aa, but at more random orbital phases for component Ac. However, for the 2012 spectrum (the top data set in Figure 9) shows Aa in an orbital phase just before the periastron and Ac close to the apastron. From a direct inspection of the overall shape and changes of the profile, we notice that while the Aa component is in an orbital phase close to the periastron, the line-peak profile remains similar in shape and amplitude. However, once the Aa component is at ${\phi }_{\mathrm{Aa}}=0.88$, the line-peak profile becomes wider. From those observational changes, we infer that most of the narrow emission, close to the line peak, comes from the Aa component.

On the other hand, changes in the blueshifted wing of Hα appear to be related to component Ac. Notice how an absorption is observed in the wing at orbital phases close to the apastron (${\phi }_{\mathrm{Ac}}=0.57$ and ${\phi }_{\mathrm{Ac}}=0.43$) and how it changes to a secondary emission peak before the periastron (${\phi }_{\mathrm{Ac}}\,=0.89$). However, besides the described changes in the line, without a detailed study of both binary components, it is still challenging to fully constrain which component of the binaries is contributing to the emission and/or absorption part of the profile. Future monitoring of the spectrum not only with GRAVITY but with other facilities like FEROS or UVES/VLT is necessary to complete the spectral analysis of the source.

In addition to the observed changes in the FEROS data, it is important to highlight that the Brγ line profile seen with GRAVITY is quite similar to the Hα line observed by Mayer et al. (2001) with the 1.4 m CAT telescope at La Silla Observatory and the bottom spectrum presented in Figure 9. Both the Brγ and the Hα lines have a peak emission of about 10%–20% above the continuum with a sharp edge at the blue part of the spectrum and an extended wing toward the red. The Brγ red wing extends up to ∼540 km s⁻¹, which is similar to the velocity of the Hα red wing reported by Mayer et al. (2001). This line feature is a very strong indication for the presence of a stellar wind component.

The X-ray luminosity measured by Parkin et al. (2011) supports the presence of shock tracers like Brγ in the wind–wind collision regions of the quadruple system. However, the semi-detached nature of the binary Ac suggests the existence of a substantial case-B mass transfer (Leung et al. 1979) in which the transferred shock-heated plasma might be the main contributor to the observed Brγ emission. One possibility for recognizing which of the two binaries is the source of the spectral features is to disentangle the individual spectra of the binaries Aa and Ac. With the interferometric observations, this is possible following the procedure detailed in Chesneau et al. (2014). This method relies, first, on the determination of the binary separation through a model fitting to all of the data (like the one we applied in Section 3.1). Second, once the separation is determined, one can analyze the relative fluxes C_i of each of the components per independent channel. Once C_i have been obtained, one can get the spectrum S_i of each component using the following expression:

$$\begin{eqnarray}&&{S}_{i}=\displaystyle \frac{{C}_{i}{S}_{* }}{{\displaystyle \sum }_{j=1}^{n}{C}_{j}},\end{eqnarray} \tag{ 5 }$$

where S_* is the spectrum of the source. After trying this method, we concluded that the data quality is not good enough to perform such an analysis with the current data sets. Two main reasons are identified for this limitation: (i) the relatively large variation of the observables in adjacent channels across the bandpass and (ii) the relative faintness of the detected lines. However, future observations with GRAVITY at the Unit Telescopes would allow us to perform this analysis.

To investigate the region in which the observed Brγ line is formed, we computed the differential phases (Millour et al. 2006). Since we observe the target as a binary, the pseudo-continuum has a contribution from both components of the binary; thus, detecting a differential signature implies a relative change in the brightness of the components at the position of the line. Nevertheless, we do not detect any systematic change in this observable at the position of Brγ compared to the pseudo-stellar continuum. As an example, Figure 10 displays the Brγ differential phases of one of our GRAVITY data sets (similar trends are observed in the other ones). Since the flux ratio of the binaries appears to be preserved across the Brγ profile, this line is produced at each of the internal binaries instead of being formed in a more extended region as a consequence of the interaction between both subsystems Aa and Ac. However, this result is constrained by the S/N in our data. Considering an average 3σ deviation of the differential phases of ${{\rm{\Delta }}}_{\phi }^{3\sigma }\sim 15^\circ$, we derived an upper limit for the line-emitting regions of p = 0.157 mas for a baseline of B = 120 m, following the next expression (Kraus 2012):

$$\begin{eqnarray}&&p=\displaystyle \frac{{{\rm{\Delta }}}_{\phi }^{3\sigma }}{2\pi }\,.\,\displaystyle \frac{\lambda }{B}.\end{eqnarray} \tag{ 6 }$$

Since the observed lines are highly variable depending on the different orbital phases of the spectroscopic binaries, future observations of this source with higher S/N are required to systematically determine the astrometric offsets associated with the different orbital phases of the inner binaries. With at least 3.5 mag more in sensitivity, GRAVITY using the 8 m Unit Telescopes will provide us with the required capabilities to (i) monitor the temporal changes in the lines, (ii) to disentangle the spectra of the spectroscopic binaries, and (iii) to provide an estimation of the orbital solution for the internal binaries (with an expected differential phase calibration of ${{\rm{\Delta }}}_{\phi }^{3\sigma }\sim 2^\circ$ or 20 μas at 120 m baselines).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This work demonstrates the enormous scientific capabilities of the new GRAVITY instrument for the characterization of multiplicity in massive systems. The system HD 93206 A is a showcase for this kind of analysis. It represents the first object of a future survey with GRAVITY. In particular, we have obtained the following results for HD 93206 A:

• 1.  
  We characterize the relative astrometry between the spectroscopic binaries Aa and Ac with a precision of tens of microarcseconds. This corresponds to a dramatic improvement of one and three orders of magnitude compared to the previous astrometric measurements with PIONIER/VLTI and NACO/SAM, respectively. The derived separation between Aa and Ac is $\rho \sim 30$ mas with a position angle P.A. ∼ 328° and ${\rm{\Delta }}K=0.32$ mag, with a detected projected motion toward the northwest from the 2012 data. This demonstrates the gravitational boundedness between the two spectroscopic binaries and provides a rough estimate of their mutual orbital period of P ≤ 60 years. The precise determination of the Aa–Ac astrometry over several future epochs (∼5 years) will be crucial to constrain the fundamental parameters of the system and its formation history, and to predict whether the system is long-term (over a few million years) stable or not.
• 2.  
  From the CANDID analysis of our GRAVITY observations, we derived a maximum magnitude difference of ${\rm{\Delta }}K\,\sim 5$ mag to detect companions, with a 3σ confidence level, in a range of 3–50 mas around the primary. This represents an increment in sensitivity of ∼2 mag from the one achieved with PIONIER/VLTI. The current surveys with the first generation of VLTI instruments have a limiting ${{\rm{\Delta }}}_{m}\sim 3\mbox{--}4$ mag between binary components with separations from 1 to 10 mas. Thus, the present result confirms GRAVITY as a unique instrument to explore a new parameter space in the characterization of massive multiple systems. Here, we have demonstrated the superb calibration and sensitivity of the single-field mode of the instrument. Nevertheless, the expected improvement in the GRAVITY performance (particularly its complete integration with the 8 m Unit Telescopes) over the next year, in combination with the so-called "dual-field" mode of the instrument, will provide us with the opportunity to look for companions in stars as faint as ∼10 mag with the ATs and ∼16 mag with the Unit Telescopes, letting us to completely characterize systems like HD 93206.
• 3.  
  A dedicated Principal Component Analysis of the GRAVITY data allowed us to extract the calibrated spectrum of the source in which shock tracers (Brγ) and stellar (He i) lines are observed. The detected spectral features have a maximum intensity ∼10% above the pseudo-stellar continuum. Complementary FIRE data supports the variability of the K-band lines, and visible FEROS spectra indicate that the daily variability is linked to the motion of the multiple components of the quadruple. It appears that changes around the peak of the Hα are associated with the Aa component, while changes in the wings correspond to the Ac component.
• 4.  
  We constrain the near-infrared wind–wind collision region associated with the observed X-ray emission of the source. From our differential phase analysis at the position of Brγ, we confirm that it is formed in very compact regions inside the inner binaries and not from a more extended region between both subsystems Aa and Ac. This is consistent with previous models that suggest that the X-ray emission is arising mostly within the spectroscopic binaries, in particular from component Ac. With the current data, we could not confirm which of the two spectroscopic binaries is the main X-ray emitter. However, from the 3σ standard deviation of the differential phases from the continuum baseline, we infer an upper limit for the line-emitting regions of p ∼ 0.157 mas (0.37 au). Future observations with the Unit Telescopes in "single-field" mode will provide us with a calibration of ${\rm{\Delta }}\phi \sim 1^\circ$ to constrain the line-emitting region up to scales of 20 μas (0.05 au at the distance of the source). If the detected changes in Hα are applicable to Brγ, detecting systematic changes in the differential phases would allow us to disentangle the astrometric motion of the individual components within the spectroscopic binaries. Monitoring projected astrometric motions combined with spectroscopic data and a proper model of the emission line variability will let us to characterize, for the first time, the fundamental parameters of the individual components in the system. This information would provide the mass distribution among the components, test coplanarity among the different orbital planes, and confront the different formation scenarios in this kind of system.

We thank the anonymous referee for useful comments to improve our manuscript. We thank the ESO and GRAVITY consortium Science Verification Team for all of their support in carrying out the presented GRAVITY observations. We thank Gabriel Ferrero for his help in processing the FIRE spectra. This research has made use of the Jean-Marie Mariotti Center SearchCal and LitPro services co-developed by FIZEAU and LAOG/IPAG, CRAL and LAGRANGE. J.S.B. acknowledges the support from the Alexander von Humboldt Foundation Fellowship programme (grant number ESP 1188300 HFST-P); this work was also partly supported by OPTICON, which is sponsored by the European Commission's FP7 Capacities programme (grant number 312430). J.M.A. acknowledges support from the Spanish Government Ministerio de Economía y Competitividad (MINECO) through grant AYA2013-40 611-P. A.A. acknowledges support from the Spanish Government Ministerio de Economía y Competitividad (MINECO) through grant AYA2015-63939-C2-1-P, co-funded with FEDER funds. F.C. acknowledges that the research leading to these results has received funding from the European Community's Seventh Framework Programme under Grant Agreement 312430 (OPTICON). The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement No. [614922].

Software: GRAVITY pipeline (Lapeyrere et al. 2014), esorex (http://www.eso.org/sci/software/cpl/esorex.html), P2VM algorithm (Tatulli et al. 2007; Lacour et al. 2008), FIREHOSE, MIDAS, emcee (Foreman-Mackey et al. 2013), Scipy (http://www.scipy.org/), LitPro (Tallon-Bosc et al. 2008), CANDID (Gallenne et al. 2015).

Figures 11 and 12 display the distributions of the V² S/N and of the errors in the closure phases after flagging