AN INTERFEROMETRIC AND SPECTROSCOPIC ANALYSIS OF THE MULTIPLE STAR SYSTEM HD 193322

Massive O-type stars are usually found with one or more nearby companions (Mason et al. 1998, 2009). Most of these luminous stars are very distant, and consequently, we generally only detect their very nearby companions through their Doppler shifts or very distant companions that are angularly resolved. We must rely on high angular resolution observations of the few nearby cases to detect those elusive, mid-range separation binary stars. One of the most revealing examples is HD 193322 (O9 V:((n)); Walborn 1972), the central star in the sparse open cluster Collinder 419. The distance to the cluster is 741 ± 36 pc according to a recent study by Roberts et al. (2010). The star's complex multiplicity became apparent with the discovery of a companion Ab through speckle interferometry observations by McAlister et al. (1987). They designated the system as CHARA 96 Aa (McAlister et al. 1989), and subsequent speckle measurements detected its orbital motion (Hartkopf et al. 1993; Hartkopf 2010). The Aa,Ab pair was also recently resolved through the technique of lucky imaging by Maíz Apellániz (2010). The composite optical spectrum is dominated by a relatively narrow-lined component Ab1, and Fullerton (1990) discovered significant radial variations in this component indicative of a spectroscopic binary. The first spectroscopic orbit for Ab1 was presented by McKibben et al. (1998), who determined an orbital period of 311 days. In addition to the close Ab1,Ab2 spectroscopic pair and the speckle Aa,Ab pair, there is another wider companion B at an angular separation of 2.68 arcsec (Turner et al. 2008). Components C and D are more distant companions that also occupy the central region of Collinder 419 (Roberts et al. 2010), and it is uncertain whether they are orbitally bound to the central multiple system. A mobile diagram presenting the known components of the system is illustrated in Figure 1. The long orbital period estimate for A,B pair is based upon the probable masses (see Table 8 below), distance, and the assumption that the projected separation is the semimajor axis.

Zoom In Zoom Out Reset image size

Measurements of the orbital motions of the stars in this system offer us the means to estimate the masses of the components. We have continued our interferometric and spectroscopic monitoring of the system over the last decade, and here we present a progress report on the orbits, mass estimates, and spectral properties of the component stars. Combined speckle and interferometric observations of the motion of the Aa,Ab pair are used in Section 2 to derive a preliminary orbit for the wide system. In Section 3, we present new long baseline interferometric measurements of the Ab1,Ab2 pair that are calibrated using the visibility of the Aa companion. In Section 4, we describe a diverse collection of spectroscopic observations which we use to derive a revised orbit for the narrow-lined Ab1 component in the close pair. In Section 5, we apply a Doppler tomography algorithm to a subset of the blue spectra to extract the spectra of the components. Finally in Section 6, we discuss the masses and other properties of the components of the system.

The orbital motion of the Aa,Ab pair has been followed since its discovery through continued speckle interferometry observations made mainly with the Mayall 4 m telescope at Kitt Peak National Observatory (KPNO; McAlister et al. 1989, 1993; Mason et al. 1998, 2009). The date, position angle θ, and separation ρ of these previously published observations are collected in Table 1 for convenience. While outliers exist, for this magnitude difference Δm and ρ regime, the errors from speckle interferometry measures are approximately 05 in position angle and 0.5% in separation. We have also measured the relative motion through optical long baseline interferometry (OLBI) with the GSU CHARA Array at Mount Wilson Observatory (ten Brummelaar et al. 2005). These separation and position angle measurements are determined by measuring the fringe packet separation when possible along two pairs of baselines with approximately orthogonal directions projected onto the sky (the separated fringe packet or SFP method; Farrington et al. 2010). This separation is determined by fringe fitting in order to avoid shifts caused by overlap. Other methods, like fits of the fringe envelope, may suffer if one fringe packet overlaps the secondary lobe of the other and causes the center of the fringe envelope to move. Like lunar occultation measurements, a single baseline measurement provides a separation in one direction only. Each of these measurements defines a line in astrometric space, and observations at several projected angles are required to fully define the position of the secondary. The location of the secondary is defined as the point with the minimum total rms distance from these lines, as weighted by the variance of the fringe separations. Formal errors are calculated using a method analogous to a χ² analysis, and the errors for ρ and θ are defined as the distance change required to increase the weighted rms by 1.0 in χ². These data have fairly low signal-to-noise ratio (S/N) and for many epochs we have data from only a single baseline with a varying position angle from diurnal motion, so there is more scatter in the resulting astrometry than in the speckle data. All the speckle and CHARA measurements are collected in Table 1.

We determined a new orbital solution for the wide pair using the combined set of speckle and long baseline interferometric observations. Note that we did not include the measurement from Maíz Apellániz (2010) using lucky imaging with the AstraLux instrument because of its relatively large error. All the measurements were initially assigned equal weight, but in the orbital fitting process we identified three discrepant points with large residuals that we subsequently zero weighted in the final fit (those dates are marked in Table 1). The orbit was determined using the grid search method described by Hartkopf et al. (1989). The orbital elements are listed in Table 2 and the appearance of the visual orbit on the sky is shown in Figure 2. The original speckle data set covers about 12% and the recent speckle and SFP data cover about 15% of the 35 year orbit. Note that we ignored making corrections for the center of light motion of the Ab binary because the largest astrometric shifts are expected to be small, ≈0.5 mas.

Zoom In Zoom Out Reset image size

The interferometric fringe patterns of the close Ab1,Ab2 pair of stars overlap even with the longest baselines available at the CHARA Array, so we cannot use the SFP method to measure the relative separation. However, the interference of the two fringe patterns of the inner pair causes a modulation of their combined visibility amplitude with changes in the projected baseline separation, and we can use this modulation to estimate the binary separation projected along the baseline position angle in the sky (Hummel et al. 1998; Boden et al. 2000; Raghavan et al. 2009). The calibration of visibility is aided when the signal of a nearby star produces an SFP that can be used to calibrate the visibility of the central binary. The details of this method are outlined by O'Brien et al. (2011).

We obtained 195 observations of HD 193322 over 24 nights between 2005 and 2010 using the CHARA Array Classic beam combiner (ten Brummelaar et al. 2005). The observations were made with the near-IR K' filter and a variety of different baselines. These measurements were a series of approximately 200 recorded fringe scans sampled at a frequency of 150 Hz. The scans were reduced by standard techniques (ten Brummelaar et al. 2005), and a subset of 100 scans with the best S/N were selected. We fit these scans with fringe patterns for both the calibrator (Aa) and target (Ab) (identified according to the orbital solution from Section 2 and the position angle of the specific baseline; see details of the procedure in O'Brien et al. 2011). We retained only those visibility measurements for which the fractional difference between a first estimate and final mean was less than 20%, and these best-case values were used to form the mean visibility for each component. Note that we used the same data set as that for the wide pair (Section 2), but due to the stronger selection criteria not all data sets yielded useful visibility amplitudes. We next determined the ratio of the mean visibilities of the target and calibrator. This observed ratio is related to the ratio of the individual visibilities for the target and calibrator by

$$\begin{equation} V_{\rm Ab} / V_{\rm Aa} = (F_{\rm Aa} / F_{\rm Ab}) V_{\rm Ab,o} / V_{\rm Aa,o}, \end{equation} \tag{ 1 }$$

where F_Aa/F_Ab is the monochromatic flux ratio in the K' band. The angular diameter of the calibrator Aa is small enough that V_Aa ≈ 1.0 for our observations, but we estimated the single-star visibilities of each component based upon the projected baseline of observation and the predicted angular diameters. A value for the flux ratio of $F_{\rm Aa} / F_{\rm Ab} = 10^{-0.4 \triangle m_{\rm wide}} = 0.92$ was adopted based upon a fit of the observed ratios and a binary model (see below), and this parameter essentially normalizes the target visibility so that the upper distribution of the visibilities has a mean of one. Note that we expected that the ratio would be F_Aa/F_Ab ⩾ 1 based on the speckle orbit assignments (Section 2), but we suspect that this difference is probably insignificant given the uncertainties in the component flux fractions. The results for the Ab pair are given in Table 3 which lists the heliocentric Julian date of observation, the corresponding orbital phase in the Ab1,Ab2 orbit (Section 4), the projected baseline B_p and position angle ψ_p of observation, the calibrated visibility and its associated error, and the observed minus calculated difference O − C in visibility from the adopted model fit. These data are available as part of the online materials in an OIFITS file.

The modulation of the visibility ratio depends on the known projected baseline length and position angle and the effective wavelength of the K' system, plus the unknown projected binary separation and magnitude differences, ▵m_close = −2.5log (F_Ab2/F_Ab1) and ▵m_wide = −2.5log (F_Aa/F_Ab). The latter magnitude difference normalizes the visibility according to the relation given above, while the former magnitude difference sets the amplitude of the visibility modulation with the baseline (Raghavan et al. 2009; O'Brien et al. 2011). Following the example of O'Brien et al. (2011), we explored the orbital parameters of the Ab1,Ab2 pair by creating a set of model visibilities for each of the observed times and baseline parameters and then forming the χ² statistic for the differences between the observed and model visibilities. The solution is found by determining the orbital parameters and magnitude differences in a high-resolution grid of values that minimize χ². For this application of the method, we set the orbital period and epoch from the spectroscopic elements for the circular orbit of Ab1,Ab2 (Section 4) and then made a grid search for the best-fit values of the angular semimajor axis a, inclination i, and longitude of the ascending node Ω, plus the magnitude differences ▵m_close and ▵m_wide.

We found that the solutions always arrived at similar estimates for the magnitude differences, so we set these magnitude differences and performed a grid search over i and a, the two parameters of physical interest. For each selection of (i, a), we determined the best fit for the sky orientation parameter Ω over the full range of values in steps of 2°. The resulting χ² estimates are plotted in a gray-scale diagram in Figure 3 in the (i, a) plane for grid increments of ▵i = 2° and ▵a = 0.05 mas. Here intensity is scaled between the lowest (black) and highest (white) χ² over the grid. If we assume a distance d from the cluster fitting results of Roberts et al. (2010), then Kepler's third law relates the known period P, the total mass, and a by

$$\begin{equation} (M({\rm Ab1})\, + \,M({\rm Ab2}))/M_\odot = {{(a d)^3} \over {P^2}} = \! \left({{a} \over {1.22 \,{\rm mas}}} \right)^3 \! \left(\!{d\over {741 \,{\rm pc}}} \! \right)^3. \end{equation} \tag{ 2 }$$

Next, we can use the spectroscopic semiamplitude K for component Ab1 (Section 4, Table 6) to derive a relation for the mass of Ab2 as a function of i and a,

$$\begin{equation} M({\rm Ab2}) = {{a^2} \over {P \sin i}} {{K (1-e^2)^{1/2}} \over {29.8 \,{\rm km \, s}^{-1}}} = 0.455 \, M_\odot {a^2 \over {\sin i}} \left(\! {d\over {741 \,{\rm pc}}} \! \right)^2. \end{equation} \tag{ 3 }$$

Then we can find the mass of Ab1 from a relation for the mass ratio,

$$\begin{eqnarray} {{M({\rm Ab1})} \over {M({\rm Ab2})}} &=& {{29.8 \,{\rm km \, s}^{-1}} \over {K}} {{a d \sin i} \over {P \sqrt{1 - e^2}}} - 1\nonumber\\ &=& {{a \sin i} \over {0.819 \,{\rm mas}}} \left({d\over {741 \,{\rm pc}}}\right) - 1. \end{eqnarray} \tag{ 4 }$$

Thus, each point in the (i, a) plane is associated with specific masses M(Ab1) and M(Ab2), and we can use the relations above to construct loci of constant primary and secondary mass in Figure 3 (shown by solid and dashed lines, respectively).

Zoom In Zoom Out Reset image size

Inspection of Figure 3 indicates that there are two broad valleys in the (i, a) plane where the fits are relatively good, one with i < 90° for counterclockwise motion in the sky and another with i > 90° for clockwise motion. Within these valleys there are three locations with comparable minima, but all of these are associated with extreme masses: (i, a) = (66°, 4.6 mas) (χ² = 210) and (i, a) = (118°, 4.7 mas) (χ² = 215) where the masses are too high and (i, a) = (38°, 3.3 mas) (χ² = 167) where the masses are too low (see Section 6 below). We think that the i < 90° valley probably represents the best family of solutions, since the trends in χ² are more or less continuous there as expected. Between a = 2.7 and 5.0 mas, the valley floor never rises above a reduced chi-square of χ²_ν = 1.18 (with 190 degrees of freedom, equal to 195 measurements minus five fitting parameters). Although a purely statistical assessment would restrict the solution space to the valley region around (i, a) = (38°, 3.3 mas), the fact that the reduced chi-square is close to unity along the length of the valley suggests that at this stage it is premature to rule out any of this solution space. In Section 6 below, we present several lines of argument that indicate that the actual solution lies in the mid-range of this valley at a = 3.85 mas (χ² = 225), so we will tentatively adopt this value and present the associated solution for the other orbital parameters in Table 2, Column 3. The errors associated with i, a, and Ω reported in Table 2 correspond to their range over the length of the valley from a = 2.7–5.0 mas. Note that because the visibility oscillation depends on the absolute value of the projected separation, there is a 180° ambiguity in our derived value of Ω. We found that the best-fit magnitude differences are ▵m_close = 2.11 ± 0.06 mag and ▵m_wide = 0.086 ± 0.012 mag (Ab brighter than Aa). In order to show how well the model and observed visibilities agree, we plot the individual and calculated visibilities for this solution for each night in Figures 4 and 5, and we find that the fits are satisfactory for most of the nights.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The projection of the orbit on the sky for this solution is illustrated in Figure 6 where filled circles indicate the calculated positions at the times of observation. The distribution of the observations in orbital phase appears to constrain the minor axis of the projected ellipse better than the major axis. With the minor axis fixed, the major axis will vary with inclination as a(minor) ∼ acos i or $a \propto \sec i$, and this relation approximately describes the position of the χ² valley in Figure 3. The orbital orientations of the wide and close orbits appear quite different (compare Figures 2 and 6), and the mutual inclination of the orbital planes ϕ is given by Fekel (1981) as

$$\begin{eqnarray} \cos \phi &=& \cos i_{\rm close} \cos i_{\rm wide}\, \pm \,\sin i_{\rm close} \sin i_{\rm wide}\nonumber\\ && \times \cos (\Omega _{\rm close}\, - \,\Omega _{\rm wide}), \end{eqnarray} \tag{ 5 }$$

where the second term may assume either sign because of the 180° ambiguity in the determination of Ω_close. The two solutions, ϕ = 38° and 85°, indicate that the two orbits are probably far from coplanar (ϕ = 0°).

Zoom In Zoom Out Reset image size

Spectroscopy potentially offers us the means to determine the masses and spectral properties of the components. However, because the stars are so close, conventional, ground-based spectroscopy records the flux of all three stars (usually plus component B at a separation of 27; Turner et al. 2008), and since their orbital Doppler shifts are comparable to the line widths, the resulting line blending problem is daunting. Nevertheless, the spectral properties of the components are sufficiently different in this case that we may attempt radial velocity measurements. The appearance of the optical spectrum is dominated by a narrow-lined component that corresponds to the primary star in the close orbit, Ab1 (McKibben et al. 1998). The secondary in the close orbit, Ab2, is fainter and contributes little to the composite spectrum (Section 5). Furthermore, there is a very broad lined component that appears to follow the motion of Aa in the wide orbit (Section 5). Because the lines of Aa are so broad and shallow, they essentially act to depress the continuum in the vicinity of the narrow lines of Ab1, and since the velocity range of Ab1 is smaller than the full width of the lines of Aa in general, the presence of the broad component has little influence on velocity measurements of Ab1 (but see a discussion of blending effects below). Here we present radial velocities for the Ab1 component and show that they represent the sum of orbital motions in both the close and wide systems.

We collected 31 new spectra for measurement from sources which are summarized in Table 4. The columns list a source number (for identification with the specific radial velocities listed in Table 5), date of observation(s), spectral range used in the measurement, the spectral resolving power, number of spectra made at that time, and the observatory, telescope aperture, and spectrograph of origin. We obtained most of these spectra in runs at the Kitt Peak National Observatory 0.9 m coudé feed and 4 m Mayall Telescopes, the 3.6 m Canada–France–Hawaii Telescope (CFHT), the 2.5 m Nordic Optical Telescope, the Lowell Observatory 1.8 m Perkins Telescope, and the Herzberg Institute of Astrophysics, Dominion Astrophysical Observatory 1.8 m Plaskett Telescope. These were augmented with publicly available spectra from the archives of the University of Toledo Ritter Observatory 1.0 m telescope (Morrison et al. 1997), the Observatoire de Haute-Provence 1.9 m telescope and ELODIE spectrograph (Moultaka et al. 2004), and the Indo-U.S. Library of Coudé Feed Stellar Spectra (Valdes et al. 2004; made with the KPNO 0.9 m coudé feed telescope). All these spectra were reduced by standard techniques and transformed to a continuum normalized flux representation on a heliocentric, log  wavelength grid. Atmospheric telluric lines were removed from the red spectra by division with a pure atmospheric spectrum. This was done by creating a library of spectra from each run of a rapidly rotating A star (usually ζ Aql), removing the broad stellar features from these, and then dividing each target spectrum by the modified atmospheric spectrum that most closely matched the target spectrum in a selected region dominated by atmospheric absorptions.

The spectra form a diverse collection with a wide range in resolving power and wavelength coverage. In order to measure radial velocities in a consistent way, we cross-correlated each of the spectra with a standard, model spectrum (rest frame) from the grid of synthetic spectra from Lanz & Hubeny (2003). From an initial inspection of the observations, we selected a model with Galactic abundances, an effective temperature of T_eff = 34.8 kK, a gravity log of g = 4.0, a projected rotational velocity of Vsin i = 50 km s⁻¹, a wavelength-dependent limb-darkening coefficient from Wade & Rucinski (1985), and an instrumental broadening appropriate for the specific observed spectrum. The cross-correlations were generally made over the wavelength range given in Table 4, although in some cases regions with strong interstellar features were omitted. The resulting cross-correlation functions were always single-peaked, and we measured the radial velocity and its associated error using the method of Zucker (2003). The results are presented in Table 5 which lists the heliocentric Julian date of mid-exposure, the corresponding orbital phases in the close and wide systems (see below), the measured radial velocity and its associated error, a correction term for line blending effects, the observed minus calculated velocity residual from the fit (see below), and the observation source number from Table 4. Note that for completeness we have included in Table 5 velocities published earlier by McKibben et al. (1998; indicated by a 0 in the final column).

Although the line blending effects from the spectral components of the other stars are generally small, they tend to bias the measurements toward the systemic velocity and lead to a slight underestimate of the orbital semiamplitude. The radial velocity offset caused by line blending will depend on the character and velocity shift of each component, the spectral features measured, and the spectral resolution of the observation. In order to make a simple correction for line blending effects we adopted the following procedure for each observation. We first determined model synthetic spectra for each stellar component (Section 5) for the spectral range and instrumental broadening of the observation. The models for components Aa, Ab2, and B were co-added according to the adopted fluxes and to the Doppler shifts for the time of observation. Then we formed a series of model spectra by adding in component Ab1 for a grid of assigned velocity offsets, and we measured the radial velocity in these composite spectra using the same cross-correlation method applied to the observations. This led to a relation between the actual and measured radial velocity for each observation, and we interpolated within this relation at the observed radial velocity to determine the offset correction for blending ΔV_r(blend) = V_r(actual) − V_r(measured), which is given in Table 5, Column 6. The average of the absolute value of the offset correction for blending is small, 2.3 km s⁻¹, but the individual offset corrections are larger for the lower resolution spectra where line blending is more severe.

The velocities of component Ab1 depend on its orbital motion in the close binary plus the motion of the Ab1,Ab2 center of mass in the orbit of the Aa,Ab system. Our first solutions for the orbital motion of the close binary clearly showed long-term variations in the residuals that followed the motion predicted for Ab in the wide orbit. Thus, we fit the observed radial velocity variations as the sum of motions in the close and wide binaries. This was done iteratively using the orbital fitting program of Morbey & Brosterhus (1974). We first made a general fit of the velocities for the close system, and then we made a constrained fit of the velocity residuals, by fixing P, e, and ω from the visual orbit of Aa,Ab to find the semiamplitude K and epoch of periastron T for the wide system. The resulting solution of the long-period orbit was then used to correct the observed velocities for motion in the wide orbit, and a new solution was found for the close orbit. This procedure quickly converged to yield the orbital elements given in Table 6. Note that we assigned each measurement a weight proportional to σ(V_r)⁻² in making the fits, and we zero weighted four measurements that had unusually large residuals from the final fit (dates indicated in Table 5 and shown as open circles in Figure 7). Table 6 lists the solutions both with and without application of the offset correction for line blending, and they are generally very similar except for the slightly larger semiamplitude K that results when accounting for line blending. Since the line blending problem is significant, we adopt the corrected velocity solutions that are illustrated in Figure 7 (close orbit) and Figure 8 (wide orbit) and that form the basis for the residuals O − C given in Table 5, Column 7. We found that the eccentricity associated with the close orbit is not statistically different from zero according to the criterion of Lucy & Sweeney (1971), so we present circular elements in Table 6 (where the epoch T is defined as the time of maximum radial velocity or, equivalently, the time of crossing the ascending node). The long orbital period of the wide system, P = 35 years, places HD 193322 among the top 1% of known spectroscopic binaries with very long periods (Pourbaix et al. 2004).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The two brightest components of HD 193322, Aa and Ab1, have very broad and very narrow spectral lines, respectively, and indeed it is these properties that can help us distinguish their different orbital motions. We show in Figure 9 CFHT spectra of the He i λ5876 profile from 1986 and 2008. During this interval, the broad component moved slightly redward as expected for the anti-phase velocity curve of Aa between wide orbit phases 0.82 and 0.44 (Figure 8). We collected all the available red spectra that recorded He i λ5876, and we formed an average spectrum for each run in order to increase the S/N of the spectra at each epoch. We then formed model spectra for each of Aa and Ab1 from the grid of Lanz & Hubeny (2003) using projected rotational velocities and model parameters optimized to match the composite profile (see below). These model profiles were fit to the observations using a nonlinear, least-squares procedure to derive the radial velocities of Aa and Ab1 at these epochs. The derived Ab1 velocities are identical within errors to the corresponding measurements given in Table 5, and the velocities for Aa are listed in Table 7. The errors associated with the velocities of Aa are large, ±20 km s⁻¹, because this component is so broad and shallow and because the shape of the red wing is sensitive to the details of the removal of the telluric features found there. Given these larger errors and the relatively small number of measurements, we made a constrained fit of the orbital radial velocity curve of Aa by setting all the parameters from the solution for Ab (Table 6) with the exception of the systemic velocity γ and semiamplitude K, and by assigning a weight to each observation proportional to the product of the spectral resolving power and the net S/N ratio in the adjoining continuum (Column 5 of Table 7). The fit (illustrated in Figure 8) yielded γ = −12 ± 8 km s⁻¹ and K = 21 ± 16 km s⁻¹, with a residual rms =21 km s⁻¹. These measurements are consistent with the expected Doppler shifts and masses for the Aa,Ab system (Section 6).

Zoom In Zoom Out Reset image size

We estimated the spectroscopic parameters for the components by first reconstructing their spectra using the orbital velocity curves and a Doppler tomography algorithm (Bagnuolo et al. 1994) and then comparing the reconstructions with models from the grids of Lanz & Hubeny (2003, 2007). We selected nine spectra from our observations that recorded the blue portion of the spectrum with a resolving power greater than 10,000 and that covered the extremes of motion in the wide and close systems (samples 10, 12, 14, 15, and 16 from Table 4). We began by running the tomography algorithm for only two components, Aa and Ab1, however, we found that the subsequent reconstructed spectrum for Aa had a composite appearance with both broad and narrow components, unlike our expectation from the He i λ5876 profiles (Figure 9). We think that this is due to flux contamination in our blue spectra from the nearby B component (B1.5 V; Vsin i ≈ 100 km s⁻¹; see McKibben et al. 1998, Figure 2; Roberts et al. 2010, Figure 1). Although component B may be a spectroscopic binary with a low semiamplitude (McKibben et al. 1998), we simply assumed that it was stationary and contributed 11% of the total flux (Roberts et al. 2010) in the next iteration of tomographic reconstruction. The power of the tomography algorithm to derive reliable and high-quality reconstructed spectra increases with the number of spectra and with the orbital velocity range and flux contribution of the components. Unfortunately, in the case of our blue spectra of HD 193322, these criteria are really only met for component Ab1. The velocity range of Aa, for example, is so small relative to its characteristic line width that the algorithm may incorrectly assign line flux between the reconstruction of Aa and the stationary component B. We dealt with this problem by starting the initial guess for components Aa, Ab2, and B with model spectra rather than assuming a flat continuum spectrum (as done for Ab1). Although the resulting solutions are guided by our assumptions, they do at least show that the observed spectra are consistent with these assumptions since otherwise the reconstructed spectra would converge to an appearance different from the initial model guesses.

We show the results of the full, four-component, tomographic reconstructions in Figure 10. These representative solutions were made using the orbital solutions from Table 6, adopting a mass ratio of M(Ab2)/M(Ab1) = 0.37 (Section 6) and the flux ratios given in Table 8. These flux ratios were calculated from the K'-band ▵m_close and ▵m_wide results (Section 3) assuming that these hot stars all contribute the same proportions to the B band. The spectroscopic parameters were determined by finding the Lanz & Hubeny (2003, 2007) model that best matched the absorption line ratios and Hγ Stark broadening (dependent on T_eff and log g) and with a Vsin i value adjusted to fit the widths of the absorption lines other than Hγ. The results are listed in Columns 2 and 3 of Table 8 for Aa and Ab1, respectively, and we estimate that the associated errors are ▵T_eff = ± 1 kK, ▵log g = ± 0.5, and ▵Vsin i = ± 40 and  ± 10 km s⁻¹ for Aa and Ab1, respectively. These parameters suggest spectral classifications of O9 Vnn and O8.5 III for Aa and Ab1, respectively, based upon the calibration of Martins et al. (2005). The "nn" suffix for the former classification indicates very broad lines. The relatively good agreement between the observed and model line depths indicates that the flux ratios from interferometry (Section 3) are fully consistent with the derived strengths of the spectroscopic features. The parameters in Table 8 for component B were taken from the work of Roberts et al. (2010), and the predicted model spectrum agrees well with the narrow, stationary spectral component from the tomographic reconstruction. The results for the faintest component, Ab2, are poorly constrained because this star contributes such a small fraction of the total flux, but its spectrum suggests an early-B, dwarf classification. We used the flux ratio between Ab2 and B and the temperature of B and the theoretical main sequence (T_eff, M_V) adopted by Roberts et al. (2010) in order to estimate the effective temperature of Ab2, assuming it is a main-sequence star. The Hγ λ4340 line is the only strong feature in the reconstructed spectrum of Ab2, and the relative weakness of the He i λλ4387, 4471 lines suggests that Ab2 may also be a rapid rotator with broad and shallow lines. Note that the Vsin i estimate for Ab2 in Table 8 is only approximate and may be subject to significant revision.

Zoom In Zoom Out Reset image size

One of the primary goals of this study was to determine the masses of the component stars. Since most of our results are preliminary, we cannot yet derive accurate masses, but the observational work does demonstrate the potential for improvement with further interferometric and spectroscopic observations. The mass sums for the wide and close systems can be determined from the angular semimajor axes and orbital periods (Table 2) and the distance d = 741 ± 36 pc for Collinder 419 (Roberts et al. 2010). The mass sums (Equation (2)) are (M(Aa) + M(Ab1) + M(Ab2))/M_☉ = (53.2 ± 11.5)(d/[741 pc])³ for the wide system and (M(Ab1) + M(Ab2))/M_☉ = 35.6(a/[4.0 mas])³(d/[741 pc])³ for the close system (where a is the angular semimajor axis of Ab1,Ab2).

To obtain the individual masses, we need to explore the solution space from the interferometric observations of the close binary (Figure 3) and from the spectroscopic orbit of Ab1 (Table 6). In particular, we can use the location of the χ² valley in the (i, a) plane of Figure 3 to derive a family of solutions based solely upon a (since i = i(a) from Figure 3). We show the derived individual masses as a function of a in Figure 11. The mass of Aa is set from the difference of the total mass of Aa,Ab and the mass of Ab1,Ab2 (from a and Equation (2)); the mass of Ab2 is from Equation (3) (from a, K(Ab1), and the i = i(a) relation in Figure 3); and the mass of Ab1 is from the difference M(Ab1, Ab2) − M(Ab2). These are all plotted in Figure 11 surrounded by a gray zone corresponding to the acceptable range in cluster distance. We see that there is a strict upper limit of a < 4.57 mas required to keep M(Aa) > 0. Furthermore, we also see that while the masses of Aa and Ab1 cover a significant range, the mass of Ab2 changes little over the range in a. This is also shown in Figure 3, where the location of the χ² valley is close to a contour of constant M(Ab2).

Zoom In Zoom Out Reset image size

Another constraint on the mass of Aa can be formed independent of the details of the Ab1,Ab2 orbit by applying Equation (3) to the wide orbit. We take a, P, e, and i from the visual orbit of the wide system (Table 2) and combine these with the orbital semiamplitude K(Ab) from spectroscopy (Table 6) to obtain M(Aa)/M_☉ = (18.7 ± 3.7)(d/[741 pc])². The  ± 1σ region from this relation is plotted as the thick line segment for M(Aa) in Figure 11, and it corresponds to a range in semimajor axis of a = 3.96 ± 0.14 mas.

It is reasonable to assume that all three stars are main-sequence objects given the position of the Aa,Ab system in the color–magnitude diagram (Roberts et al. 2010). The K-band fluxes of massive main-sequence stars scale with mass M as F∝M^2.30 for stars in this mass range according to the models of Marigo et al. (2008), so we can use this relation to predict the flux ratio between any pair of stars according to the mass relations shown in Figure 11. The positions where the model flux ratios match the observed ones (Table 8) are indicated by pairs of symbols in Figure 11. These flux ratio relations indicate a semimajor axis range of a = 3.80 ± 0.08 mas.

A final constraint can be set from the overall fluxes and absolute magnitude of the combined system. ten Brummelaar et al. (2000) estimate that the apparent V-band magnitude of Aa,Ab is V = 5.96 ± 0.02 mag, and using the distance and extinction for the star from Roberts et al. (2010), we estimate that the absolute magnitude is M_V = −4.35 ± 0.12 mag. The individual absolute magnitudes of the components are given in Table 8. We can apply the mass–absolute magnitude relation (M, M_V) for the main sequence from the models of Marigo et al. (2008) to obtain a prediction for the total absolute magnitude for each of the mass combinations shown in Figure 11. The best match occurs for the masses obtained at a = 3.81 mas, approximately where Aa and Ab1 have equal masses. For those masses the predicted absolute magnitude is M_V = −5.12, which is significantly brighter than the estimate above from observations. The models predict even brighter fluxes if either Aa or Ab1 is more massive, and the horizontal line at the bottom of Figure 11 shows the range in a where the combined magnitude is within 0.1 mag of the faint limit.

The average estimate for the semimajor axis from the above three constraints is a = 3.85 ± 0.09 mas, and we adopt the associated mass solution as best representing the current observational data. The masses and other properties summarized in Table 8 are generally in agreement with expectations for hot, main-sequence stars (Martins et al. 2005). However, a number of significant discrepancies remain that deserve further investigation. The mass of component Ab1 is similar to that expected for an O8.5 III star (≈24 M_☉; Martins et al. 2005), but the star's absolute magnitude is about 1.8 mag fainter than typical for such stars. This discrepancy hints that Ab1 may be a dwarf rather than a giant star. The overall faintness of the system compared with expectations for the stars' masses may indicate that the distance estimate needs to be revised downward (leading to lower masses) and/or the extinction estimate revised upward. Some confusion also remains about whether Aa or Ab1,Ab2 is brighter. As we noted in Section 3, the visibility analysis indicates that ▵m_wide = m(Aa) − m(Ab) = 0.086 ± 0.012 mag (Ab brighter than Aa), which agrees within errors with high angular resolution measurements by Maíz Apellániz (2010) using the AstraLux camera, m(Aa) − m(Ab) = 0.04 ± 0.19 mag. Although these results indicate that Ab is somewhat brighter than Aa, we refrain from re-designating the identities of the components to avoid confusion with published results.

We expect that these lingering problems will be resolved with future interferometric observations that will better sample the orbital phases of the close pair near the nodal crossings and will lead to improved constraints on the angular semimajor axis a. In addition, we clearly need to continue the long-term high-resolution work on the wide orbit to cover the missing orbital phases (Figure 2). We plan to obtain these measurements through continuing observations of this system using the CHARA Array interferometer. Additional high S/N and high-resolution spectroscopy holds the promise of delivering better orbital constraints on Aa and Ab2 that would then allow us to estimate the masses without relying on the distance of the cluster. Indeed, reliable orbital elements would render it possible to set an independent estimate of the cluster distance. We anticipate expanding the spectroscopic coverage over the next decade.

We close with some speculative remarks about the angular momentum distribution of the stars of HD 193322. It is remarkable that this system contains both a very rapidly rotating star (Aa) and a very slowly rotating star (Ab1). It is possible that Ab1 has a rotational axis with a low inclination, so that its equatorial rotational velocity is close to typical values. However, the very large line broadening of Aa places it among the most rapidly rotating O-type stars known (Penny 1996). It is possible that the angular momentum of the natal cloud led directly to rapid rotation in the case of Aa and to the formation of a binary in the case of Ab. Alternatively, there may have been some very close gravitational encounters in the early life of the system. In some circumstances, a close encounter between a binary and a third interloper can lead to a merger of two of the components and ejection of the third (Gaburov et al. 2010). It is possible that the rapid rotator Aa is such a merger product and that the runaway star 68 Cygni was the object ejected from the system (Schilbach & Röser 2008). If so, then the orbital and spin properties of the stars of HD 193322 offer key evidence about the early dynamical processes in this cluster.

The CHARA Array, operated by Georgia State University, was built with funding provided by the National Science Foundation, Georgia State University, the W. M. Keck Foundation, and the David and Lucile Packard Foundation. This material is based upon work supported by the National Science Foundation under grants AST-0606861, AST-0606958, AST-0908253, and AST-1009080. Institutional support has been provided from the GSU College of Arts and Sciences and from the Research Program Enhancement fund of the Board of Regents of the University System of Georgia, administered through the GSU Office of the Vice President for Research. B.D.M. and W.I.H. have been supported by the National Aeronautics and Space Administration under reimbursable no. NNH06AD70I, issued through the Terrestrial Planet Finder Foundation Science program. Thanks are also extended to Ken Johnston and the U. S. Naval Observatory for their continued support of the Double Star Program. M.V.M. thanks Lehigh University for an institutional grant. A portion of the research in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We are grateful to Dr. Gregg Wade and the MiMeS consortium for sharing with us their spectral data on HD 193322 in advance of publication. This work is partially based on spectral data retrieved from the ELODIE archive at Observatoire de Haute-Provence (OHP). Additional spectroscopic data were retrieved from Ritter Observatory's public archive, which was supported by the National Science Foundation Program for Research and Education with Small Telescopes (NSF-PREST) under grant AST-0440784.