IMAGING THE ALGOL TRIPLE SYSTEM IN THE H BAND WITH THE CHARA INTERFEROMETER

Algol (β Per, HD 19356) is one of the most extensively studied stellar systems in the history of astronomy. It is a triple hierarchical system where an inner eclipsing binary is orbited by a more distant outer star. The inner system of orbital period 2.87 days is a prototypical semi-detached binary in nearly circular orbit. The primary, Algol A, is a main-sequence star of type B8V. The secondary, Algol B, is a subgiant of type K2IV. As Algol A and Algol B are thought to have formed roughly at the same time, one would expect the more massive star of the pair to become a subgiant first. However, this supposition is contradicted by all observations, showing Algol B to be lighter than Algol A. The solution to this so called "Algol paradox" is that most of the original mass from Algol B, which fills its Roche lobe, is thought to have been accreted by Algol A. As the inner pair is tidally locked, this creates a strong magnetic dynamo effect and non-thermal emission observed in X-ray and radio. Using radio interferometry, Lestrade et al. (1993) detected positional displacement during the orbital revolution of the inner pair and identified Algol B as the source of radio emission. In X-rays, giant flares have been detected on Algol B (White et al. 1986; Schmitt & Favata 1999; Schmitt et al. 2003). Their sizes were estimated to a few tenths of the Algol B radius, and they were shown to follow active and quiescent periodic cycles of 49 days (Richards et al. 2003). Hα profiles provide evidence for mass transfer through a stream of gas from the Roche lobe of Algol B to Algol A (Gillet et al. 1989). In parallel, observations with EXOSAT (Favata et al. 2000), XMM (Schmitt & Ness 2004), and Chandra (Drake 2003) have shown evidence of an extended corona around Algol B. The presence of a large, permanent coronal loop around Algol B, oriented toward Algol A, has very recently been confirmed by radio images from the HSA (Peterson et al. 2010).

The tertiary component, Algol C, was first discovered by radial velocity measurement (Curtiss 1908) and revolves around the inner pair over a period of about 680 days without eclipsing. It is an Am star with metal and hydrogen lines presenting F1V characteristics (Richards 1993). Earlier studies based on astrometric, photometric, and radial velocity data allowed the determination of the masses of the inner stars (Hill et al. 1971; Tomkin & Lambert 1978) and of Algol C (Bachmann & Hershey 1975). Then light curve meta-analyses combining previous visible, infrared, and ultraviolet observations (Soderhjelm 1980; Richards & Mochnacki 1988) produced estimates of the orbital elements for the inner and outer orbits, as well as for the surface temperatures, log g, radii, and masses of the three stars. With the advent of speckle interferometry, it became possible to spatially resolve the outer component (Labeyrie et al. 1974; McAlister & Degioia 1979), leading to the determination of orbital parameters for the outer orbit (Bonneau 1979). Optical interferometry then improved upon these results with the Mark III interferometer (Pan et al. 1993). However, due to 180° ambiguities in early speckle observations, all studies based on speckle results (Bonneau 1979; Soderhjelm 1980; Pan et al. 1993) produced ascending node estimates for the inner or outer orbits in disagreement with the polarimetric (Rudy & Kemp 1978) and radio (Lestrade et al. 1993) observations. Therefore, as for given inclinations, the angle between the inner and outer orbits solely depends on the difference of the ascending nodes; this ambiguity implied that the orbits were found either roughly perpendicular or coplanar. As coplanarity was ruled out by the shallow eclipse depth of photometric observations (Soderhjelm 1980), ad hoc corrections had to be made to the ascending nodes of speckle-based studies. More recently, Csizmadia et al. (2009) resolved the inner pair with long baseline interferometry at the Center for High Angular Resolution Astronomy (CHARA), using the CLASSIC instrument working in the near-infrared K_s band, in combination with very long baseline interferometry radio observations. They determined the orbit of the inner binary to be in prograde rotation, in disagreement with the retrograde movement found previously in radio by Lestrade et al. (1993). Zavala et al. (2010) recently solved these long-standing issues, by simultaneously resolving all three stellar components in the optical with the Navy Optical Interferometer (NOI). Using angular referencing, the most precise inner and outer orbits to date were derived, with positional errors on the stellar components of the order of a milliarcsecond. The outer orbit was shown to be prograde and the inner orbit retrograde. The latest Very Long Baseline Array radio measurements are in agreement with this orientation (Peterson et al. 2011). Zavala et al. (2010) also presented the first reconstructed interferometric image of all three stars, though the short baseline lengths used at NOI precluded full separation of the inner pair components.

In this paper, we present the results of our Algol observations with the Michigan InfraRed Combiner (MIRC) in the H band, making use of the long baselines available at CHARA. In Section 3, the relative positions of the three components are determined by model fitting for the 55 epochs of our data set. New estimates of the orbital parameters are then derived in Section 4 using these stellar positions. Finally, images of the inner orbit are reconstructed for each epoch in Section 5, where we discuss the potential evidence for the phenomena detected by radio and X-ray instruments.

Our observations of Algol were conducted between 2006 and 2010 at the Georgia State University CHARA Array interferometer using the MIRC instrument. The CHARA Array, located on Mount Wilson, CA, consists of six 1 m telescopes and is the largest optical/IR interferometer array in the world (ten Brummelaar et al. 2005). Its 15 baselines range from 34 m to 331 m, providing resolutions up to ∼0.5 mas at H band and ∼0.7 mas at K band. Thus, the data collected constitute the highest angular resolution data available to date on Algol. MIRC was used to combine four CHARA telescopes together for true interferometric imaging in H band, providing six visibilities, four closure phases and four triple amplitudes simultaneously in eight narrow spectral bands for each snapshot (Monnier et al. 2006).

Our observations total 20 nights, making the accumulated interferometric data on Algol one of the largest data sets to date at CHARA. In particular, the data cover most phases of the inner pair eclipse. A complete observing log is listed in Table 1, including our calibrators. The spacing of Algol observations is not regular: Algol was sometimes observed only as a backup target, when atmospheric conditions did not allow for other fainter targets. The duration of most of these observations generally did not exceed 1 hr. Two array configurations were used, S1-E1-W1-W2 and S2-E2-W1-W2. The full complement of fringes was always obtained except on 2009 August 18 when bad seeing prevented the acquisition of the E1 fringes. Both configurations are well adapted to imaging, with roughly equal Fourier coverage in all directions. Figure 1 presents the typical Fourier coverage achieved within one of our final data sets.

Zoom In Zoom Out Reset image size

The standard observing procedures were followed, and the data were reduced using the standard MIRC pipeline described in Monnier et al. (2007). After frame co-adding, background subtraction, and a Fourier transform of the raw data, fringe amplitudes and phases are used to form squared visibilities and triple products. On all data prior to 2010, the photometric calibrations are estimated using shutter matrix measurements and partial beam chopping. Starting in 2010 data, the calibration relies on the newly available (and far superior) photometric channel technique, allowing flux measurement in parallel with the fringes (Che et al. 2010). During the next step, calibrators are used to calibrate the drifts in overall system response. Table 2 presents the sizes of the calibrators, as estimated from our ongoing study of "good" interferometric calibrators and based on three independent photometry methods, adopting different color–magnitude relations to compute the stellar sizes. It has to be noted here that 37 And is currently suspected to be a binary, though with a flux ratio between components large enough not to significantly affect the calibration. After calibration, the pipeline outputs the final calibrated power spectra and bispectra into files compliant with the OIFITS standard (Pauls et al. 2005).

Based on the orbital parameters obtained by previous studies, we expect significant stellar movement on short timescales. The orbital period of the inner pair is about 2.87 days, meaning that Algol B is moving relatively to Algol A by about 0.07 mas in an hour. The movement of Algol C also impacts the visibilities through high order lobes, which oscillate rapidly with its separation from the inner pair (ranging from 15 to 70 mas), though this phenomenon is somewhat mitigated by its much longer orbital period (about 680 days).

While model fitting can possibly be made to account for these time dependencies, image reconstruction packages currently cannot. At best, an image reconstructed using our data accumulated over an hour would be blurred, and at worse it would not converge at all toward a meaningful solution. Thus, the data have to be split into shorter sequences of no more than 10 minutes, finding a compromise between reasonable Fourier coverage and temporal consistency. This splitting is applied to about 60% of the nights and results in a total of 55 different files, each of which contains about 240 power spectra and 160 bispectra, sufficient for imaging without blurring. Figures 2 and 3 illustrate how the data accumulated on 2009 August 12 have been split into chunks.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The purpose of model fitting here is to derive the position, shape, and flux parameters of all three Algol components. These parameters will then be used to determine the orbital parameters of the inner orbit (noted A–B) and outer orbit (noted AB–C), as well as to generate realistic prior images for later image reconstruction. Our general procedure is to minimize the residuals between the interferometric observables derived from our current model and the corresponding interferometric data, using a Levenberg–Marquardt algorithm. Two possible paths were initially examined: either directly fit the orbital parameters on the full data set or estimate a separate set of parameters for every single epoch and then fit the orbital parameters using these positions. As Algol B is non-spherical (because it is thought to fill its Roche lobe), the two-dimensional projection of its surface onto the observing plane varies from night to night (Richards 1992). The "global fit" approach would require a complex three-dimensional model of the three stars, while the "night-to-night fit" approach we finally settled for is more flexible, as well as more manageable in terms of computing power and memory requirements.

Each star is modeled here as an ellipse of uniform brightness distribution. As the angular sizes of the components are only few times the array resolution, the use of limb-darkening coefficients would have a nearly negligible impact on the visibilities. In fact, limb-darkened models are not strictly necessary for orbit fitting and even detrimental as they increase the number of parameters to estimate. Furthermore we expect image reconstruction to retrieve the true flux distribution. A major effect to take into account when modeling the system is the bandwidth smearing effect on Algol C, i.e., Algol C's contribution to the visibilities is affected by its large separation from the delay-tracking center (Algol A) and the width of each of the MIRC spectral channels. With this effect, the total number of parameters to estimate for each epoch is 23: the principal axes and angles for Algol B and Algol C relative to Algol A (the available interferometric data does not give access to absolute positions), the relative flux contributions of Algol A and Algol B to the total flux (normalized to unity), the shape parameters of the three elliptical disks (semimajor axis, ellipticity, and angular orientation), and the bandwidth smearing factors for Algol C (one for each spectral channel).

The χ² formed from power spectra and bispectra is known to be multimodal: As the Levenberg–Marquardt is a gradient algorithm, the parameters have to be initialized close to their true value to prevent falling into local minima. Attempting to estimate the 23 unknowns simultaneously would only lead to inconsistent results. We determined that the positions of Algol B and Algol C relative to Algol A constitute the parameters with the most influence on the χ². Therefore our solution to this convergence problem involves simultaneous systematic grid searches on both Algol B and C Cartesian coordinates (relative to Algol A), coupled to Levenberg–Marquardt optimizations on the remaining parameters. In addition, some parameters can safely be fixed, based on the literature of the expected physical properties of the system. As Algol A and Algol B are thought to be tidally locked (Richards 1992), the inner orbit is posed circular. The principal axis of Algol B is supposed to be along the direction A–B, as expected for a Roche-lobe-filling secondary. Algol A and Algol C are both assumed perfectly circular, because they are not rapid rotators and are not filling their Roche lobes. With these assumptions the χ² fit on the remaining number of parameters is then well constrained as each data set contains between 120 and 240 power spectra as well as 60 and 160 bispectra (triple amplitudes and closure phases). Out of the eclipse periods, lower and higher bounds are imposed on each parameter solely to enforce the physical meaning of parameters. During the eclipses, the aspect ratios of Algol B are kept small (<1.05).

The fitting is done in two steps—first a coarse grid search on Algol B and Algol C positions, then a finer one. The coarse search uses a 10 × 10 grid with spacing 0.1 mas aimed at determining the approximate positions of the components. For each epoch, the center of the Algol C grid is initialized at its expected position based on the orbital elements from the NOI results (Zavala et al. 2010), while the center of the Algol B grid search is simply set to the origin (0, 0) (i.e., Algol A position). For a given night, this setup of the coarse grid search covers all the probable positions of the components. Then the finer grid search refines the estimates, using a 10 × 10 grid with spacing 0.02 mas around the most likely positions found during the coarse search. To insure convergence during both the coarse and fine searches, the fluxes and angular diameters of the star are always optimized first, while the secondary ellipticity and the bandwidth smearing factors are kept fixed to their initial estimates; then the constraints on these are relaxed and all parameters are re-optimized.

To compute the visibility contribution of Algol C, an analytic expression of the complex visibilities of a bandwidth-smeared ellipse is directly employed. For the visibility contribution of the inner pair, the presence of epochs where one component eclipses the other prevents the use of analytic expressions. Algol A and Algol B are therefore modeled as ellipses on an image array of 64 by 64 pixels with a pixel size of 0.05 mas. To prevent problems linked to the discretization of the ellipse sizes with this pixellation, the edges of the ellipses are smoothed by a sigmoid function, so that the χ² varies continuously with the size parameters. The complex visibilities are then computed by applying a Discrete Fourier Transform to the image. Finally the model power spectra and bispectra are derived from the visibilities, and the reduced χ² is computed.

Table 3 reports these reduced χ² values for every night. Typically the agreement between the model and the data is good, especially considering the unrealistic assumption of uniform brightnesses: The values lie mostly within the 1.0–6.0 range, with higher χ² values encountered during the eclipse epochs where this model breaks down. Figure 4 presents the residuals of a typical fit, which show no noticeable trend.

Zoom In Zoom Out Reset image size

The full output of our model-fitting software consists in reduced χ² maps as a function of Algol B and Algol C positions, as presented in Figure 5. The most probable positions of Algol B and Algol C correspond to the peak value of these maps. For each map, the contour at 1σ (68% confidence level) is computed and fitted by an ellipse. The positions of the peaks for Algol B and Algol C as well as the parameters of the error ellipses (major and minor axes, angles) are recorded in Table 3.

Zoom In Zoom Out Reset image size

The fits also produce estimates for the angular diameters of all three stars. The estimates are ϕ_A = 0.88 ± 0.05 mas, ϕ_B = 1.12 ± 0.07 mas, and ϕ_C = 0.56 ± 0.10 mas. For Algol B, this angular diameter corresponds to the major axis of the ellipse (the aspect ratio of the Algol B ellipse changes with the phase of the inner orbit, between 1.04 and 1.22 for our data sets, which we attribute the projection of the Roche lobe).

Concerning Algol C, we note that the difficulty of estimating the bandwidth smearing factors is reflected here as a larger error bar associated with its angular diameter. Our angular diameters are not directly comparable with the published values in Table 1 of Zavala et al. (2010), as these suffered from a conversion error (i.e., the parallax values in milliarsceconds and parsecs have been confused). Nevertheless, as they were based on the linear diameters estimated in Richards (1993), we settle for comparison with Richards (1993) only. We will assume the most recent value published for the orbital parallax, 34.7 ± 0.6 mas, independently determined by optical (Zavala et al. 2010) and radio (Peterson et al. 2011) interferometry, and more precise than the latest Hipparcos estimate of 36.17 ± 1.40 mas (van Leeuwen 2009). Finally the linear stellar radii computed with this parallax are given in Table 4: Our results are consistent with the radius estimates derived by light curve analysis in Richards (1993), though they admittedly have greater error bars due to the night-by-night nature of our fits and the uniform brightness assumption.

Using the estimates of the relative positions of the components established in Section 3, the orbital parameters for the inner and outer orbits are derived.

The orbital elements for both orbits and the mass ratio m_A/m_B are fitted simultaneously, again using the Levenberg–Marquardt algorithm. The parameters are initialized around the published values in Zavala et al. (2010), for which the corresponding reduced χ² is about 13.0. After convergence our final best χ² is 0.94. In addition to these orbital parameters, Kepler's third law applied to both the inner and outer orbits allows us to derive the mass ratio m_C/(m_A + m_B) independently from parallax assumptions. The parameters corresponding to the best fit are given in Tables 4–6, with error bars corresponding to a 68% confidence level. The probability distributions of the parameters were obtained with the bootstrapping-with-replacement technique (using 50, 000 bootstrap samples) and thus takes into account correlated errors. In the case where distributions were found to be roughly symmetrical about the peak values, only a single error bar is given. Figures 6 and 7 show, respectively, the inner and outer orbits corresponding to our best fit, as well as the estimated positions of Algol B and Algol C as determined in the previous section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most of our estimated orbital elements are in agreement with previous results based on light curve analysis (Richards 1993) and radio and optical interferometry (Zavala et al. 2010; Peterson et al. 2011). In particular, the respective prograde and retrograde movements of the outer and inner orbits are unambiguously verified (see also the reconstructed image sequence of the inner orbit in the next section). The mutual inclination angle i_m between both orbits is given by

$$\begin{equation} \cos i_m = \cos i_1 \cos i_2 + \sin i_1 \sin i_2 \cos (\Omega _1 - \Omega _2), \end{equation} \tag{ 1 }$$

where i₁ and i₂ are the inclination of the inner and outer orbits, and Ω₁ and Ω₂ are the associated ascending node angles. The mutual inclination determines the period of orbital precession between the orbits, and consequently the amplitude of the variations of inclination of the close binary due to dynamical effects (Borkovits et al. 2004). Previously published values close to 100° (Kiseleva et al. 1998) would imply a detectable change of inclination of the close binary from 1° to 3° in the last century, which in turn would lead to observable variations in the minimum depth of the light curves of the eclipse and the disappearance of these eclipses after a few centuries (Csizmadia et al. 2009). This is however contrary to the actual observations: Algol eclipses have been known since Antiquity and therefore the minimum depth is not thought to have significantly changed within error bars. The mutual inclination estimate derived from our bootstrap results is i_m = 9020 ± 032. It is completely consistent with photometric observations, and much closer to exact perpendicularity than the most recent estimate of 95° ± 5° (Zavala et al. 2010).

The main point of disagreement between previous studies and our work is that we found the semimajor axis of the inner orbit to be slightly smaller, in turn implying a smaller value for the total mass of the inner binary. However, to disregard this shorter axis result we would have to invoke a 15% error in our estimation of the position of the center of mass of the secondary. As our analysis is essentially based on the determination of photocenters, a potential source of error to consider is the asymmetric brightness distribution on the surface of a Roche-lobe-filling star. We quantified the magnitude of such an effect, by simulating the equatorial brightness profile of Algol B, taking into account both gravity darkening and limb darkening as shown in Figure 8. The positions of the photocenter and of the center of mass were found to coincide to 1%. Arguably the existence of an additional strong proximity effect (the primary heating of the atmosphere of the secondary) could also explain a photocenter shift. Detailed modeling of such an effect, however, goes beyond the scope of this paper, requiring complementary data (imaging in Section 5 based on current data did not allow detection of the proximity effect). Finally, the semimajor axis of the inner binary was also found to be 15% smaller than previous optical results by recent radio (Peterson et al. 2011) and X-ray observations (Chung et al. 2004), a result strikingly similar to ours. While both of these papers hypothesized the existence of a coronal component around Algol B, no coronal signature has been detected in the near-infrared yet. Our result thus suggests the possibility that the center of light and center of mass actually roughly coincide, with previous optical results suffering from imprecisions. Our observations indeed benefited from much higher angular resolution than previous interferometric studies and from the recent disambiguation of the orbital angles from Zavala et al. (2010).

Zoom In Zoom Out Reset image size

Consequently, we recommend that our orbits (like all previously published orbits based on interferometric data) should be interpreted as photocenter orbits. We also note that we did not notice significant differences between the centers of light estimated by model-fitting uniform ellipses and estimated by using the real brightness distribution. In Section 5, we demonstrate that the semimajor axes estimated from the disk model and from the reconstructed images differ by no more than 5%.

The inner binary mass ratio m_A/m_B = 4.56 ± 0.34 agrees with the most recent radial velocity measurements (Hill et al. 1993, 1971; Tomkin & Lambert 1978), which measured m_A/m_B = 4.6 ± 0.1. Our larger error bar is due to the lesser intrinsic precision of interferometry in this parameter, as the mass ratio influences the position of Algol C only by a small amount. Figure 9 shows that the χ² of the Algol C orbit fit possesses a very flat minimum valley, resulting in the large error on the mass ratio.

Zoom In Zoom Out Reset image size

With the assumption of a parallax of 34.7 ± 0.6 mas (Zavala et al. 2010; Peterson et al. 2011), the masses of the three stars are completely determined and we report the results in Table 4. Algol A and Algol B are lighter than that found by Zavala et al. (2010) and Richards (1993). We note here that the tertiary mass published in Zavala et al. (2010) is numerically inconsistent with the corresponding semi-axes and periods, which seems to be due to slightly different parallax values adopted for the inner and outer orbit fits (C. Hummel 2011, private communication). With our hierarchical fit, we find the mass of Algol C to be larger than that published in Richards (1993). Overall our mass estimates are consistent with all measurements (radial velocity and system dynamics) and give a sensibly different view of the system. In particular the picture of the mass transfer within the inner binary may have to be modified.

The amount of data present in each MIRC data sets corresponds to a dozen minutes of observation each night. All data sets except one (2009 August 18) contain four-telescope data. Consequently there is enough phase information to attempt "model-independent image reconstruction" on each night, where the prefix "model-independent" underlines that the method does not rely on a specific astrophysical model.

The image reconstruction procedure belongs to the class of "ill-posed" inverse problems, as it attempts to reconstruct an image conventionally consisting in a large number of pixels (typically a few thousand) using the less numerous interferometric data points (typically a few hundred). Expressing the problem in a Bayesian framework shows that the solution can be given by the regularized maximum likelihood method (Baron et al. 2010; Thiébaut & Giovannelli 2010). The target image is the array of fluxes $\widehat{\boldsymbol{x}}=\lbrace x_0, \ldots x_{n-1}\rbrace$ that minimizes the sum of a term linked to the data (the likelihood) and a term reflecting all prior information (the regularization) under the two constraints of image positivity and of normalization of the image to unity, i.e., in the mathematical form

$$\begin{equation} \widehat{\boldsymbol{x}} = \underset{\boldsymbol{x} \in \mathbb {R}^n}{{\rm arg\,min}} \lbrace \chi ^2 (\boldsymbol{x}) + \mu R(\boldsymbol{x}) \rbrace \end{equation} \tag{ 2 }$$

subject to positivity (∀i, x_i ⩾ 0) and normalization (∑^{n − 1}_{i = 0}x_i = 1). The likelihood term measures the distance between the observed power spectrum and bispectrum and the same quantities derived from the tentative image. Minimizing this term enforces the presence of flux as indicated by the data, but the minimization presents many local minima. The addition of the regularization function $R(\boldsymbol{x})$ discriminates between these local minima, as well as prevents overfitting, detrimental to image quality. The factor μ in Equation (2) controls the relative weight of the χ² and regularization terms, and can be chosen so that the actual reduced χ² is roughly unity for the reconstructed image (the selection of the truly optimal μ is a difficult problem which goes beyond the scope of this paper). Most regularizers are computed in the image plane, but also possess good frequency extrapolation properties in the Fourier plane. Hence they allow reconstructions to routinely achieve super-resolution, i.e., to produce images with an effective resolution typically about four times greater than the physical array resolution.

The image reconstruction software SQUEEZE (Baron et al. 2010) has been used to obtain the reconstruction presented in this paper. It combines the Markov Chain Monte Carlo and gradient-descent approaches used, respectively, by its predecessors MACIM (Ireland 2006) and BSMEM (Baron & Young 2008). SQUEEZE offers a vast choice of regularizers, ranging from classic entropy to novel wavelet regularizations. As the images of Algol are expected to consist of three compact limb-darkened stars, we choose to use the total variation regularizer, an edge-preserving regularizer which ensures both the smoothness and the compactness of the stars while strongly penalizing stray flux (Rudin et al. 1992). In our implementation, the total variation is implemented as

$$\begin{equation} R_1(\boldsymbol{x}) = \sum _i - {\rm TV}(x_i) = \sum _i |x_{i+1}-x_i|. \end{equation} \tag{ 3 }$$

A second regularizer is also employed in the form of a prior image to initialize the Markov Chain to a sensible starting point as well as to prevent the exploration of pixels where the presence of flux is known to be very improbable. The prior image $\boldsymbol{m}$ is created by convolving the best model image of Section 3 to one-fourth of the array resolution. The entropic prior expression is the Burg entropy (Burg 1975):

$$\begin{equation} R_2(\boldsymbol{x}) = \sum _i - \log (x_i/m_i). \end{equation} \tag{ 4 }$$

The number of pixel elements in the Markov Chain is set to 4000, and the length of the Chain to 10, 000 iterations. A common size and pixellation is chosen for all 55 images to allow easier comparison. The highest instrumental resolution is given by the largest CHARA baseline (S1-W1, 330 m), corresponding to 0.5 mas in the H band. Taking into account the expected amount of super-resolution, our choice of pixellation for the reconstructed images is thus 0.1 mas. Meanwhile the distance from Algol C to the inner pair is at most 100 mas, thus requiring an image size of 1000 × 1000 pixels. This is a factor 100 more pixels than typical binary reconstructions from non-simulated data in optical/infrared interferometry (Zhao et al. 2008). Nevertheless, all 55 reconstructions converged without notable issues, with final reduced χ² between 1.0 and 2.0.

Figure 10 presents the reconstruction for 2009 August 12. The three stars are present within a small field of 16 mas, Algol C being close to its periapsis, and they are well resolved.

Zoom In Zoom Out Reset image size

Algol A appears as a bright and near-perfect circular disk with a small amount of limb darkening. Algol B is fainter and elongated toward Algol A, as expected due to it filling its Roche lobe. Algol C also appears elongated in the direction of Algol A–B, though this is only here a purely optical-numerical artifact due to bandwidth smearing.

All the reconstructed images of the inner pair are presented in Figure 11, with epochs sorted by increasing phase of the inner eclipse. An animation generated from this sequence is also available in the online version of this paper. These images should be interpreted with caution, due to the comparatively poor resolution of the interferometer compared to the size of the stars, and to the typically limited dynamic range of image reconstruction with four telescopes (about 10:1). Overall the aspect of Algol A barely changes from image to image, while Algol B's elongation varies as expected with the phase of the eclipse.

Zoom In Zoom Out Reset image size

Due to the sporadic nature of our Algol observations, the phase coverage is unequal and incomplete, and in particular there are no data available beyond phase 0.861. About 20% of the images correspond to snapshots taken during the eclipse epochs. Images with a phase below 0.08 correspond to the primary eclipse, when Algol A is occulted by Algol B. The progressive separation of the stars is clear as the phase increases up to 0.25 (even though the data for these images were taken several years apart). Images with a phase between 0.45 and 0.55 correspond to the secondary eclipse. During part of the secondary eclipse, Algol B seems to completely disappears behind Algol A despite its larger size, an effect we attribute to its lower brightness and the low dynamic range.

In a significant number of images the secondary is clearly limb brightened, a phenomenon we attribute to gravity darkening. We can neither rule out nor confirm the presence of spots on the stellar surfaces, such as those hypothesized in Richards (1990, 1992). Similarly, we are unable to confirm the existence of a stream of matter between the stars or of an active corona around Algol B such as detected in radio (while the 2009 August 18 image at phase 0.463 would seem to show some interesting interaction between the stellar pair, it is the least reliable of all images due to data availability on three telescopes only).

Using our reconstructed images, we can now quantify the impact of our assumption of uniformly bright stars adopted in Section 3. As a function of the inner orbit phase, Figure 12 presents the amount of correction to the principal axis of A–B introduced by using the brightness distributions given by the imaging. No significant trend is detected. The maximum correction is 0.057 mas, roughly corresponding to half a pixel, and thus always within one standard deviation of the positions given in Table 3. Consequently, this validates our modeling in Section 3 and increase our confidence in the derived orbital solution.

Zoom In Zoom Out Reset image size

From the analysis of the CHARA/MIRC infrared data on Algol we derived new orbital elements of the triple system with unprecedented precision. The respectively prograde and retrograde nature of the outer and inner orbits is unambiguously confirmed and the mutual inclination of the orbital planes is found to be extremely close to perfect perpendicularity. Our results suggest that the semimajor axis of the inner orbit is shorter than previously reported by other methods. While this effect could in theory be due to a shift of the photocenter with respect to the center of mass, we did not find evidence to support this assumption. If our orbits reflect the center of mass orbit, then the masses of the primary and secondary are found to be lighter than suggested by previous works. Using model-independent image reconstruction techniques, we reconstructed a sequence of 55 images of the inner binary with an effective resolution of 0.2 mas. This "movie" of Algol A and Algol B covers most phases of the inner pair, including the primary and secondary eclipses. While the array resolution did not allow us to confirm the detection of potential flares or spots on the surface of the stars, the Roche-lobe shape of the secondary is clearly visible.

To further improve the orbital solution presented in this paper, three-dimensional and time-dependent image reconstruction algorithms are currently being developed at the University of Michigan. By modeling the secondary as a Roche lobe and by imaging directly on spheroids, we expect to considerably enhance the quality of the reconstructions, enabling us to detect potential surface features on the stars.

We thank Bob Zavala (NOI) and Christian Hummel (ESO) for their constructive criticisms on the manuscript. We are grateful to the reviewers for directing our attention to the potential effects of the asymmetric brightness distribution on a Roche-lobe surface. The authors acknowledge funding from the NSF through awards AST-0807577 to the University of Michigan. Operational support for the CHARA Array is provided by the National Science Foundation through grant AST-0908253 and by the College of Arts and Sciences at Georgia State University.

Facility: CHARA (MIRC) - Center for High Angular Resolution Astronomy interferometer