HUBBLE SPACE TELESCOPE ASTROMETRY OF THE PROCYON SYSTEM^*

Figure 1 illustrates a typical pairing of a short unsaturated exposure in F218W for astrometry of Procyon A, and a much longer exposure at the same pointing used to analyze Procyon B. The inset in the center shows Procyon A from a 0.11-s exposure, superposed on a 100-s frame in which the WD is easily detected in spite of the neighboring, grossly overexposed image of the primary star.

Zoom In Zoom Out Reset image size

In our first visit in 1997, we were too aggressive in choosing an exposure time for A of 0.14 s, resulting in its image containing saturated pixels in most of the exposures. Fortunately, due to dithering, there were two short-exposure frames in which A remained unsaturated. Even at the reduced 0.12 s used in the 1998 visit, two of the short exposures were again saturated for Procyon A. For the remainder of the WFPC2 observations, we set the short exposures to the WFPC2 minimum of 0.11 s, and none of them were saturated.

In all cases we used the individual c0m.fits images from the archive pipeline for the astrometric analysis. These frames have bias subtraction and flat-fielding applied, but do not include any cosmic-ray removal, geometric correction, or drizzle processing. Each short- and long-exposure pair was taken at a different dither position, using fractional-pixel offsets to sample the point-spread function (PSF), plus shifts of a few integer pixels to average out the impact of detector defects (such as hot pixels). We checked for discrepant measurements due to cosmic-ray impacts, but found no cases where they had caused a problem in our relatively short integrations.

In the analysis of the 1995 WFPC2 observation by G00, the relative positions of Procyon A and B were determined by cross-correlation of the short-exposure image of A with the long-exposure image of B. The uncertainty for this measurement was at a level of about ±0.2 pixels (∼0009). However, the accumulation of WFPC2 data from several programs, including ours, which used F218W between 1994 and 2009, makes it possible to apply a more precise astrometric analysis based on PSF fitting. (We did try the original cross-correlation approach for our WFPC2 data, but found that the errors were about 50% greater than those based on PSF fitting.)

PSF fitting is based on an empirically derived, over-sampled representation of the image structure, obtained by combining numerous high-S/N exposures taken at many independent pointings. As indicated in Table 1, we have 53 unsaturated exposures of 0.11–0.14 s on Procyon A. In addition, we included 96 archival observations of the standard star Grw +70°5824, obtained for WFPC2 F218W photometric calibrations between 1994 and 2009, and two frames of the standard star BD +17°4708 from 2004. That made a total of 151 well-exposed F218W images of stars having colors similar to Procyon A and B, positioned near the center of the PC chip, for input to the PSF determination.

The approach we used is described in Equations (2) and (3) of Gilliland et al. (1999). A uniform spatial grid on a scale finer than the native pixel size is first defined. (In this case, we chose a factor of 50 finer than the input scale, because it allowed re-use of existing codes developed for Kepler analyses, but the results are insensitive to the exact choice.)

The individual pixel values, after normalization of all inputs to unit volume, are accumulated into a weighted sum at each over-sampled grid point, using a Gaussian weighting based on separations of each input pixel from the accumulation grid. The width of the Gaussian weighting function is a free parameter. Adoption of too small a Gaussian width results in over-fitting the data, while too wide a weighting function suppresses available resolution in the resulting PSF. We used a Gaussian weighting width of 0.416 pixels full width half maximum (FWHM), which minimized the scatter in measurements at the same epoch.

Developing the over-sampled PSF requires precise knowledge of the relative centering of each input image. The solution is therefore iterative, since precise relative positions are best determined through fitting the over-sampled PSF to individual images. Fortunately a simple first-moment estimate of image positions for all inputs is accurate enough to start a rapidly convergent iterative cycle of determining an over-sampled PSF, revising the image positions, and recalculating the PSF.

With the PSF defined, we then obtained the relative positions of individual images of Procyon A and B by fitting a bi-cubic interpolation function (Press et al. 1992) to the PSF, and then employed a nonlinear least-squares fit (Bevington 1969) to determine the relative x, y centers of both stars. These fits used the central 21 pixels (5 × 5 box without corners), after experiments showed that smaller or larger fit domains performed marginally less well.

Procyon B lies in the extended wings of the PSF of Procyon A. Simulations indicated that these wings shift the measured position of B by less than 0001 at times of greatest separation early in the WFPC2 series, increasing to 0003–0004 near closest approach. We therefore derived a deep PSF by stacking all of the strongly saturated Procyon A images, which we then subtracted before performing PSF fits for the position of B.

WFPC2 had significant geometric distortion, due both to the camera optics and a manufacturing defect in the CCDs. We applied geometric-correction terms for the optical distortions, and the "34th-row" detector defect, from Gonzaga & Biretta (2010). The geometrically corrected x, y positions were then converted to angular units using an F218W plate scale of 0045437 pixel⁻¹, with a nominal fractional error of ±0.0003, adjusted slightly for differential velocity aberration (using the VAFACTOR keyword in the image headers), all as described by Gonzaga & Biretta.¹²

Because of our technique of short exposures followed by long exposures, our measurements of the separation and position angle are subject to a systematic offset due to telescope pointing drift (which occurs even when the telescope is locked on guide stars). Gilliland (2005) showed that drifts of 0010–0015 are typical during HST orbits, which translates to about 00025 for a pair of A, B exposures taken over about one-fifth of an HST visibility period.

Since our images only show the two components of Procyon, we are unable to establish firmly that the pointing drift was always present, and if so what its direction was; but the effect of such drifts is likely to dominate over other terms (such as influence of differing stellar colors or changing telescope focus on the PSF, plate-scale changes due to telescope "breathing," residual contamination from component A, etc.). We have therefore estimated the errors of the average positions at each epoch by combining in quadrature the standard error based on the observed measurement scatter with a systematic term of ±00025 for telescope drift. (Although drift is systematic within a single HST visit, a range of different telescope orientations was used across the different epochs, so it is appropriate to treat drift error as a random term.)

Lastly, we determined the absolute J2000 position angle of B relative to A, using the ORIENTAT keyword in the image headers, which gives the orientation on the sky of the image y axis. The error on position angle includes two terms. The first arises from the errors of derived x, y positions of A and B, estimated as described above. This term in position angle will be inversely proportional to the lever arm provided by the changing separation of A and B during their orbit. A second term arises from uncertainties in the absolute HST roll angle. We assume a 1'' error on guide-star positions, observed with the FGS over a ∼1000'' baseline, which translates to an angular error of ±0028. We combine these two uncertainties in quadrature. The final astrometric results from the WFPC2 F218W images are given in Table 3, lines 1 and 4 through 14.

The HST archive contains images of Procyon in the WFPC2 F1042M filter obtained at two epochs in 1997. The primary aim of these observations was a search for faint companions of 23 nearby stars, including Procyon, but no new companions were found by the proposing team (Schroeder et al. 2000). For our astrometric measurements, these frames raised the challenge that Procyon B is nicely exposed, but the core of the image of A is strongly saturated. This forced us to develop a means of using the outer portions of the PSF of Procyon A, where the spatial information is dominated by the diffraction spikes, to obtain its centroid location.

We faced a similar problem in our complementary HST program on the Sirius system, to be discussed in a separate forthcoming paper. Our WFPC2 images of Sirius were taken in the same F1042M filter, and are likewise saturated for Sirius A. In order to test methods for centroiding saturated images, we carried out a calibration program (Program ID: CAL/WFPC2-11509) on the star 109 Virginis (spectral type A0 V, V = 3.73). This star has a color similar to that of Sirius, and not extremely different from Procyon. It is sufficiently faint that unsaturated images in F1042M can be obtained in short exposures (0.23 s), along with saturated images from longer integrations (600 s). We obtained a set of three dithered pairs of short and long exposures on this star.¹³

We initially considered an approach for astrometric analysis of the saturated images of Procyon in which we would develop an over-sampled PSF for the unsaturated regions of the images, including especially the diffraction spikes. As in the case of the unsaturated F218W images described above, development of the PSF for saturated F1042M images requires an ensemble of input data. Table 4 lists the images of Sirius and 109 Vir that we used (in addition to the two sets of Procyon A F1042M images listed in Table 1).

However, we found that the appearance of the diffraction spikes is unstable. Figure 2 shows examples of the variable spike structures in saturated images of Sirius and Procyon in the F1042M filter. The intensities of the spikes vary by large amounts as functions of distance from the center of the stellar image in a quasi-periodic fashion, which does not reproduce well from epoch to epoch. The structure of these intensity variations appears to depend strongly on small differences in the location of the star in the field of view. Therefore we did not see any straightforward means of defining an over-sampled PSF for the unsaturated outer regions of the deep exposures.

Zoom In Zoom Out Reset image size

We instead adopted an alternative approach of fitting straight lines to the diffraction spikes, and determining the image centroid from their intersection point. Our procedure was to estimate the location of the pixel nearest the center of the saturated image, and then search inwards toward this point along each of the four diffraction spikes until the first saturated pixel was encountered. From that pixel outward, we calculated the sums of intensities along each diffraction-spike axis and along the two neighboring parallel axes one pixel away on either side. These sums were accumulated in two sequential segments, each 30 pixels long, for a total length of 60 pixels. Then a parabola was fit to the three sums; the peak of this parabola marked the location of the diffraction-spike axis in the direction orthogonal to the spike. We found good consistency between results from the first and second 30-pixel segments along the spikes, and thus combined them to form a single center for each of the four spikes. The intersection point of the lines connecting the symmetric diffraction-spike centers then defined the stellar centroid.

Application of this approach to the three pairs of calibration frames on 109 Vir gave mean offsets between the spike intersection point and the PSF-derived centroid of −0.004 ± 0.002 pixels in x, and +0.044 ± 0.040 pixels in y (or, in arcseconds, −00002 ± 00001 and +00020 ± 00018, respectively; the errors are estimated from the scatter among the three measurements). We simply adopted these as (small) corrections to be added to the A-component positions from the long exposures.

To determine the centroids of Procyon B from the images, we followed the approach used for the F218W data. That is, an over-sampled PSF was first derived, using the same Gaussian weighting approach. We used 69 individual inputs, listed in Tables 1 and 4, consisting of all deep exposures of Procyon B and Sirius B in F1042M, plus the three unsaturated exposures on 109 Vir. Background removal was done by taking a median of an annulus from 13 to 23 pixels out from B, then subtracting it to remove the pedestal due to the wings of A.

Due to a MgF₂ lens immediately in front of the CCDs, there is a weak dependence of the WFPC2 plate scale on wavelength. A plate scale for F1042M images is not provided in Gonzaga & Biretta (2010). However, by plotting the plate scales listed by Gonzaga & Biretta against the index of refraction of MgF₂ at the effective wavelength of each filter, we found a tight, linear correlation. Only a slight extrapolation to the wavelength of F1042M was needed to estimate its plate scale. We adopted a relative plate scale of 1.00048 compared to the fiducial F555W value, for a net of 0045577 pixel⁻¹ for the PC chip.

Typical WFPC2 F1042M visits consisted of five or six exposures of about 8 s, and another five or six of about 60 s. Thus there are two distinct clumps of medium and long integrations. To see if there was a dependence of the astrometric results on exposure time, we compiled means and scatters within the medium and long blocks separately for each epoch. The average difference in measured separations from the medium- and long-exposure sets was an inconsequential 00007. The average scatter of separations was 00026 within the medium exposures, and 00035 for the long exposures. We concluded that we could safely combine the results from all of the exposures within each epoch.

The remainder of the F1042M astrometric analysis proceeded as described above for the F218W images. However, the Procyon astrometry is not directly affected by HST pointing drift, since A and B are measured on the same frames. On the other hand, our results hinge on the single-orbit calibration using three pairs of short and long 109 Vir exposures, which were subject to the drift error. The canonical drift allowance of 0010–0015 per HST orbit visibility translates in this case into a potential systematic error of ∼0004. We applied this value in quadrature with the the standard error based on random scatter. The astrometric results from the WFPC2 F1042M images are given in Table 3, lines 2 and 3.

The WFC3 observations of Procyon in F953N are similar to those in WFPC2 F1042M: Procyon A is saturated in all images. We adopted a similar approach for the analysis, beginning by assembling a set of images for PSF determination and a study of the use of diffraction spikes for centroiding. In addition to the WFC3 F953N images of Procyon, listed in Table 2, we have observed Sirius in this WFC3 filter in our complementary program on that binary. And we likewise carried out calibration observations (Program ID: GO-12598), in which we obtained both saturated and unsaturated WFC3 frames in F953N of the Pleiades main-sequence star HD 23886 (spectral type A3 V, V = 8.01). Table 5 lists the images of Sirius and HD 23886 that we used for these studies.

The data for all three targets were acquired using four-point dithering with the WFC3-UVIS-DITHER-BOX pattern. In most cases repeats were used at each setting within the pattern, providing 32 total exposures during each HST visit for Procyon, 28 for Sirius, and 8 for HD 23886.

The HD 23886 calibration observations, and most of our Procyon observations, used Chip 2 of the WFC3 camera, with UVIS2-C512C-SUB, the 512 × 512-pixel subarray nearest the Amp C readout. This subarray has been shown (Gilliland et al. 2010) to be the best behaved for photometry near and beyond saturation. The Sirius observations all used the larger UVIS2-C1K1C-SUB 1024 × 1024 subarray, also in Chip 2. However, our first WFC3 visit for Procyon in 2010 was obtained using the corresponding Chip 1 subarray. Lacking any supporting calibration observations for this chip, we have omitted the 2010 data from our analysis.

For the WFC3 astrometric analysis, we use the default drizzle-combined drz.fits images from the archive pipeline. These frames are created by combining the individual dithered exposures, and are fully processed to bias-subtracted, flat-fielded, and geometrically corrected images with cosmic rays removed.

The two left-hand panels in Figure 3 show two representative images of Procyon from our WFC3 F953N observations. In contrast to the WFPC2 frames in F1042M, the diffraction spikes have a smooth fall-off in intensity with radius, without any quasi-periodic fluctuations. Moreover, the image structure appears to be consistent over all of the epochs.

Zoom In Zoom Out Reset image size

We therefore developed a deep PSF, using all of the F953N exposures in Tables 2 and 5, except for the short unsaturated exposure on HD 23886. This PSF extends out to a large enough radius always to cover the location of Procyon B; a small region around B was set to zero weight in each individual exposure contributing to the deep PSF. For the weighting we adopted a Gaussian weighting FWHM of 0.832 pixels. With the deep PSF determined, we then subtracted it before using data on B (both for development of the unsaturated, core PSF, and for the subsequent centroiding of B).

The two right-hand panels in Figure 3 show the result of subtracting the best-fit deep PSF from two individual images. A region around the primary-star charge bleeds is not handled well, but this is inconsequential. Apart from this, within the diffraction spikes, and generally for fine structures in the PSF from A, the subtraction effectively removes 85%–90% of the flux.

As shown by the red boxes in Figure 3, we used only relatively small regions containing high-S/N but unsaturated point-like structures, in each of the four diffraction spikes, for the PSF fitting. The signal within these boxes exceeds that from B by over an order of magnitude.

Positions of B in all cases were determined using the PSF-fitting approach adopted for F218W. To create the PSF for this purpose, we started by stacking the 25 drizzled images of Procyon B, Sirius B, and HD 23886 (unsaturated). For the Procyon B and Sirius B inputs, the underlying light from A was first subtracted, using the deep over-sampled PSF of A.

Since we used two independent methods for fitting A and B in the same frames, and the PSFs do not have absolute centroids, it is important to apply a calibration using the images of HD 23886. We found that corrections of +0.468 pixels in x, and +0.449 pixels in y, needed to be added to the A-centroid technique results to bring them into alignment with the B technique. This leaves the possibility of telescope drift during the calibration observations unaccounted for. Fortunately, however, inspection of the HD 23886 images showed that it has a (previously unknown) faint companion, offset by ∼20 pixels in x, and ∼5 in y, from the bright star, which is detected in both the short and long exposures.¹⁴ This allows a direct correction for HST pointing drift (particularly valuable in this case because one of the short/long pairs was split across two spacecraft orbits). Monitoring the faint-star position with unsaturated PSF fits in both the long and short exposures indicated a correction for drift of +0.026 and −0.018 pixels in x and y, respectively, in the above sense. Thus the net calibration zero-point corrections are +0.494 and +0.431 pixels in x and y. Since the companion is faint in the short exposures, the precision is rather low, and we adopt a systematic error term of ±0004 for the WFC3 astrometry within epochs. This was added in quadrature as a random term since visits are at effectively random orientations.

Although the pipeline images are geometrically corrected, at the time of our initial analyses the geometric distortion in F953N had not been calibrated as well as for the more frequently used WFC3 filters. In particular, it was not included in the study of WFC3 plate scales by Kozhurina-Platais (2014). However, a search of the HST archive yielded a set of frames in F953N of the cluster ω Centauri. We performed a new analysis of these images, generating new geometric-distortion calibration reference files paralleling those in the work just cited for other filters. After these files were incorporated into the calibration database at STScI, we retrieved the data again, and the results presented here make use of the new calibrations. The final plate scale adopted in the pipeline reductions for F953N is 003962 pixel⁻¹. The results from this method are shown in Table 3 as the "WFC3 F953N PSF fit" entries.

Having developed the alternate technique of fitting straight lines to the diffraction spikes, and then taking the intersection of these lines as the centroid of component A for the F1042M data, we also applied this technique to the WFC3 F953N data. We again used the short and long exposures on HD 23886 to calibrate the offset between the spike-determined position of A and the PSF-determined position of B. In this case, the corrections are +0.055 and +0.102 pixels, to be applied to the position of A. The results from this method are labeled in Table 3 as "WFC3 F953N spike fit" values. The largest absolute difference in A–B separation between the two methods over the four epochs is 5 mas, with a mean of 1.6 mas, and a standard deviation of 2.9 mas. This suggests that the two techniques yield comparable results, with differences between them consistent with the stated error bars. We therefore averaged the results, and show them at the bottom of Table 3, labeled "WFC3 F953N average."

Our HST measurements of the Procyon system are extremely precise, compared to ground-based data, but they cover less than half of only one orbital period. Thus the historical ground-based data are important in constraining the orbital elements, especially the period. The available visual observations of Procyon, from 1896 to 1932, were assembled by S51 (his Table 8). Since 1932, according to the Washington Double Star (WDS) Catalog maintained at the USNO, there have been only nine further published measurements of Procyon. Three of these measurements are from HST observations in 1995 and 1997 (G00; Schroeder et al. 2000), now superseded by our present results, and leaving only six new published ground-based observations since 1932.

The early observers would often report measurements averaged over several observations taken over relatively short intervals. Occasionally, the observer would recompute the averages in a subsequent publication based on a different combination of the measurements. This sometimes led to redundant listings for the same measurements. We cross-compared the S51 tabulation and WDS catalog with the original publications, and adopted the values published most recently by the observers. Additionally, in compiling his data, S51 sometimes averaged observations that had not been averaged by the original observers, and that are now listed individually in the WDS. In these instances, we adopted the individual measurements as listed in the WDS. Table 6 indicates the observations that we removed from the S51 and WDS listings, and which measurements we used to replace them.

Table 7 gives the complete list of edited ground-based measurements that were initially used in our orbit fit. Some of the observations were badly discrepant and were removed in our final fit; these are identified in the table by a superscript c in the first column. Our fitting procedure and rejection process are described below in Section 4.2. In his tabulation, S51 had corrected the position angles for precession to the J2000 equinox; in our Table 7 we have similarly corrected the position angles for the ground-based measurements after 1932 to the J2000 equinox (except for CCD and adaptive-optics observations, which we assumed to be reported for J2000).

We fitted a visual orbit simultaneously to the HST and ground-based measurements (Tables 3 and 7 respectively; for the HST WFC3 data, we used the "F953N average" values). We used a Newton–Raphson method to minimize χ² by calculating a first-order Taylor expansion for the equations of orbital motion. For the HST data we used the measurement errors directly from Table 3 in computing χ². The ground-based observers typically did not estimate errors for their measurements, so we adopted an iterative approach to optimize the weighting of the ground-based data in our orbit fit, and to reject outliers. In the first step of the iterative procedure, we fit an orbit to the ground-based data only and applied uniform uncertainties to these measurements to force the reduced ${\chi }_{\nu }^{2}$ to equal unity (where ν is the degree of freedom). In the second step, we used these scaled uncertainties to fit an orbit simultaneously to the ground-based and HST measurements. We used a sigma-clipping algorithm to reject any ground-based data point whose residual was more than three times the standard deviation of the residuals for the full data set. We repeated this procedure until no additional data points were rejected. The final data set contained the 57 measurements listed in Tables 3 and 7 (18 from HST and 39 ground-based retained in the solution). The adopted uncertainties for the ground-based separations were ±0187; we propagated this value to the position angle by assuming equal uncertainties in the R.A. and decl. directions. The historical measurements removed from the fit through sigma clipping are flagged in Table 7. Many of the rejected observations were made between 1914 and 1929, when the visual measurements were extremely difficult—or even, as suggested by S51, of doubtful reality.

Table 8 lists the final parameters for the visual orbit. The uncertainties were computed from the diagonal elements of the covariance matrix. We also investigated a solution using only the HST measurements; this solution produced uncertainties averaging about 60%–70% larger than those presented in Table 8, with the error in the orbital period more than doubled. An additional, and probably final, HST observation will be scheduled in 2016, but we expect that the historical ground-based data will continue to be an essential part of the best orbital solution.

In Figure 4 we plot the data points, both HST and ground-based, and the orbital fit. The positions of Procyon B that are predicted from our orbital elements in Table 8 are marked with open blue circles, for the HST observations only. At the scale of Figure 4, the observed HST data (filled black circles) are so precise that they appear to lie exactly at the centers of the open blue circles. For a better visualization of the errors, the two panels of Figure 5 show the residuals of the HST observations from the positions predicted by our orbital elements, in R.A. and decl. The units are now milliarcseconds (mas), rather than the arcseconds of Figure 4. The error bars are those given in Table 3, converted from separation and position angle to R.A. and decl. Based on the residual plots, there is no evidence within those errors for perturbations of the orbit by a third body. (We return to this point in Section 6.)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the elements of the relative orbit listed in Table 8, we need two further quantities in order to determine dynamical masses for both stars: the absolute parallax of the system, and the semimajor axis of the absolute motion of Procyon A on the sky.

There are three recent independent determinations of the parallax: (1) G00 obtained it from measurements of a series of ∼50 plates taken at the USNO 1.55m reflector between 1985 and 1990; (2) the parallax was measured by Hipparcos (we use the value from the new reduction by van Leeuwen 2007); and (3) it was measured with the Allegheny MAP by Gatewood & Han (2006). These results are in good agreement, and we adopt a weighted mean of the parallax values, as given in the top part of Table 9.

G00 determined the semimajor axis of Procyon A's motion from ∼600 exposures on plates obtained at six different observatories, from 1912 to 1995. We have adopted their result, a_A = 1232 ± 0008, as given in the bottom of Table 9. It agrees fairly well with a value of 1217 ± 0003 obtained by S51 from a subset of the same photographic material.¹⁵ The G00 result differs by a larger amount from the 1179 ± 0011 found by I92 from a combined analysis of RVs and a re-analysis of S51's astrometry.

Our decision to adopt the G00 value of a_A over that of S51—despite his quoted uncertainty being smaller—is based on the significant advantages of the G00 study. These include the use of plates spanning twice the time baseline (∼80 years versus ∼40 years); digital centering with a laser-encoded PDS microdensitometer as opposed to visual centering with a single-screw measuring engine; and computer-calculated plate transformations using scores of reference stars compared to the four reference stars used in Strand's manual calculations. It is possible that this last limitation of having to rely on just four reference stars might have caused Strand to underestimate the uncertainty in a_A. As a check on the uncertainty estimated by G00, Elliott Horch kindly reprocessed the 593 measures from G00 using his independent orbit-element code. The uncertainty in a_A was confirmed to be ±0008.

There is also the question of the sensitivity of the value of a_A derived by G00 to the adopted orbital elements, given our new and more precise determination of those elements. We investigated this by assuming the values in our Table 8 for all elements except the semimajor axis, and then reprocessing the photographic measures of G00, solving only for a_A. The result is unchanged and robust, with values of a_A ranging from 1231 to 1234, depending on the degree of "outlier" trimming. For these reasons, we have adopted the values of a_A and its uncertainty as given by G00.

We did not use RV information in our orbital solution, which was based purely on astrometric data. The RVs, however, provide a useful check on the validity of our final results.

Our orbital elements, along with the parallax and the semimajor axis of Procyon A's motion, allow us to predict the RV of Procyon A, apart from a constant offset due to the center-of-mass motion of the binary system. In Figure 6 (top), we compare our predictions with absolute RVs published by I92 (from photographic spectrograms, 1909–1985), and a single absolute RV by Mosser et al. (2008, from RV speedometer). Also plotted are relative RVs by Innis et al. (1994; from RV spectrometer data, 1986–1990). Innis et al. adjusted their velocity zero-point so that their RVs would match the I92 orbit predictions in the mean. We have applied I92's gamma-velocity of −4.115 km s⁻¹ to our predicted RVs.

Zoom In Zoom Out Reset image size

I92 also published a separate set of precise RVs of Procyon A, measured using a hydrogen-fluoride absorption cell, obtained over the interval 1980–1991. These velocities are on a relative scale. In Figure 6 (bottom), we compare the RVs predicted by our orbital elements with these velocities; we have arbitrarily shifted the zero-point of our predictions to match the measurements in the mean. Both of these figures show that our parameters of the Procyon system are able to predict the RV measurements very well.

Table 10 lists the dynamical masses that result from our adopted parameters. We used the usual formulae for the total system mass, $M={M}_{A}+{M}_{B}={a}^{3}/({\pi }^{3}\;{P}^{2}),$ and for the individual masses, ${M}_{A}=M(1-{a}_{A}/a)$ and ${M}_{B}=M\;{a}_{A}/a\;;$ in these equations the masses are in M_⊙, a and π in arcseconds, and P in years.

In Table 11 we present the error budgets for the masses of Procyon A and B, based on the adopted random uncertainties of each of the parameters.¹⁶ For Procyon A, the mass error is dominated almost entirely by the uncertainty in the parallax. For the WD companion, Procyon B, the parallax is again responsible for the majority of the error, but the uncertainty in the semimajor axis of Procyon A's astrometric motion also contributes significantly. Unfortunately, the mass uncertainties are unlikely to be reduced in the near future, because Procyon is too bright for its parallax to be measured by the Gaia mission (D. Pourbaix 2015, private communication).¹⁷

With the advent of asteroseismology, Procyon A was recognized as a unique object for exploring non-radial stellar oscillations. Oscillation frequencies are particularly sensitive to boundaries between radiative and convective regions. Stars near the main sequence in the mass range near that of Procyon A are believed to exhibit a convective core, a radiative envelope, and a very thin outer convection zone (e.g., Guenther & Demarque 1993).

Helioseismology shows that diffusion of helium and heavy elements in the solar interior significantly affects the solar oscillation frequencies, and it is expected also to play a role in the radiative envelope of Procyon A. In addition, the convective core overshoot at the core's edge must be taken into account. The amount of core overshoot in a star of this mass is not well known. It strongly affects the morphology of the evolutionary track and the evolutionary rate in the post-main-sequence phase of evolution where Procyon lies.

The oscillation spectrum of Procyon A has been obtained from ground-based RV measurements (Arentoft et al. 2008; Bedding et al. 2010), and from intensity observations with the space mission MOST (see Guenther et al. 2008). A Bayesian statistical study of the asteroseismic data, based on a large grid of stellar-evolution tracks, was carried out by Guenther et al. (2014; hereafter GDG14). Their tracks spanned the mass range 1.41 to 1.55 M_⊙. Other grid parameters (see Table 2 of GDG14) covered the following variables: (1) the helium and heavy-element contents by mass ranged from Y = 0.26 to 0.31 and Z = 0.014 to 0.031; (2) the mixing-length parameter, α, in the thin outer convection zone ranged from 1.7 to 2.5; and (3) the core-overshoot parameter, β, initially ranged from 0.0 to 1.0 times the local pressure scale height, H_p. Because of unanticipated evidence for large convective overshoot, the grid was eventually extended to β values as large as 2.0 H_p. Three quantities were selected as priors in the calculations, namely the mass of Procyon A (from G00), and its position ($\mathrm{log}L/{L}_{\odot }$ and $\mathrm{log}{T}_{{\rm{eff}}}$) in the HR diagram (HRD; see Table 1 of GDG14).

The strongest result of the GDG14 analysis was that all of the most probable theoretical models (with or without core overshoot, with adiabatic or non-adiabatic model frequencies, with or without diffusion in the radiative envelope, and including or not including priors for the observed HRD location and mass) were found to have masses within 1σ of the mass inferred from observations of 1.497 ± 0.037 M_⊙, as published by G00. The error bar for the most probable theoretical mass is still large, and more precise oscillation frequencies will be needed in the future to reduce it. But it is encouraging that this result is in such good agreement with the dynamical mass of 1.478 ± 0.012 M_⊙ derived from our observational analysis in the present paper.

Another result of the Bayesian analysis is relevant. The most probable models were characterized by substantial overmixing beyond the formal boundary of the convective core, with values of β as high as 1.0 H_p or even larger. This result exceeds the expected value of β = 0.2 or less, generally accepted for core convective overshoot in similar stars. This may be evidence for diffusive mixing beyond the standard overshoot region, as recently discussed by Moravveji et al. (2015) in the case of the more massive star KIC 10526294. A full understanding of this result will require continued seismic monitoring of Procyon A to improve the precision of the oscillation frequencies, as well as more sophisticated modeling in the overshoot region.

Knowing the age of Procyon A is critical to understanding the past evolution of the binary system. In conjunction with the cooling age of the companion WD, the age of Procyon A allows us to estimate the original mass of the Procyon B progenitor (see L13 and the discussion in Section 8.2 below).

We constructed grids of stellar evolutionary tracks for stars with the dynamical mass of 1.478 M_⊙ derived in the present paper, following them from the zero-age main sequence (ZAMS) to the subgiant branch. The tracks were calculated under the assumption of single-star evolution, i.e., no interaction between Procyon A and its companion during the course of its evolution from the ZAMS to the present.

We calculated models that either include or ignore the effects of element diffusion, and for amounts of convective-overshoot efficiency at the edge of the convective core of β = 0, 0.2, and 1.0 H_p. A standard model of convective core overshoot was adopted, as described in GDG14. The temperature gradient in the overshoot region was constrained to be the local radiative temperature gradient. This situation has been described as "overmixing," as opposed to "penetrative convection," where the temperature gradient is adiabatic (see Zahn 1991; GDG14). We used a near-solar metallicity of Z = 0.02, and adjusted the hydrogen abundance in the ZAMS starting model so as to ensure that each track passed through the Procyon A error box (from Doğan et al. 2010) in the HRD.

Figure 8 displays three of these tracks, constructed with the stellar-evolution code YREC (Demarque et al. 2008). All of these models include helium and heavy-element diffusion in the radiative envelope, using the formalism of Bahcall & Loeb (1990). Core overshoot is the major uncertainty in determining the ages of these models. The track plotted in red assumes no core overshoot (β = 0); the track in green has a standard value of core overshoot (β = 0.2); and the track in blue has large core overshoot (β = 1.0). The hydrogen contents in the ZAMS starting models were X = 0.672 (red track), 0.680 (green track), and 0.716 (blue track). These abundances are all consistent with accepted uncertainties in X.

Zoom In Zoom Out Reset image size

The red and green tracks are morphologically similar to tracks used in previous studies of Procyon A. In particular, note Procyon A's position just below the well-known leftward "hook," due to core hydrogen exhaustion. This near coincidence was, until the advent of precision asteroseismology, a major source of ambiguity in identifying the precise evolutionary status of Procyon A. However, seismology clearly placed Procyon A in the core-burning phase of evolution; but as discussed above it also surprisingly revealed the presence of extensive mixing in the interior outside the convective core (see GDG14). Due to the very large amount of overshoot, the blue track has a quite different morphology from the other two, and it also evolves more slowly. In this case, Procyon A lies well before the hydrogen-exhaustion phase.

The Procyon A ages based on the red, green, and blue tracks are 1.673, 1.817, and 2.703 Gyr. If one accepts the main results from the GDG14 Bayesian statistical analysis of the seismic data, then the preferred age is close to 2.70 Gyr. Such an age may be an upper limit. While a minimum age near 1.8 Gyr (as found by both L13 and GDG14) seems well established, there remains an uncertainty in the maximum age, which depends sensitively not only on the amount of chemical mixing from the core but also on the composition profile and structure above the core edge in the envelope. The recently published seismological study by Moravveji et al. (2015) of KIC 10526294, a 3 M_⊙ star near the main sequence, shows that the frequencies observed by the Kepler mission can be tightly fitted to a diffusion model in the overshoot region. Improved oscillation frequencies of the same quality will be needed to produce a similar result for Procyon A. Finally, it should also be emphasized that the validity of this discussion rests upon the assumption of single-star evolution at constant mass for the A component (see below).

Procyon B lies mostly hidden in the glare of its much brighter primary star. Virtually all that was known about it in the pre-HST era was based only on its astrometric properties and approximate brightness estimates. However, using HST, Provencal et al. (1997) acquired WFPC2 images of Procyon B in several wide- and narrow-band filters, with effective wavelengths ranging from 1600 to 7828 Å, and covering most of the star's spectral energy distribution (SED). From these photometric data they deduced a helium-composition photosphere, and estimated the star's effective temperature (${T}_{{\rm{eff}}}$ = 8688 ± 200 K) and radius (R_B = 0.0096 ± 0.0005 R_⊙). Based on the I92 astrometric mass of 0.622 ± 0.023 M_⊙, this relatively small radius called into question the assumption of the CO degenerate core that would be expected for a WD of this mass. Provencal et al. instead suggested the remarkable possibility of an iron core—placing the star in an "iron box"—as implied by the zero-temperature WD mass–radius relations of Hamada & Salpeter (1961).

The nature of Procyon B became clearer five years later when P02 used HST to obtain a series of Space Telescope Imaging Spectrograph (STIS) spectra, covering 1800–10,000 Å. These revealed the presence of C₂ Swan bands, as well as absorption features due to C i, Mg ii, Ca ii, and Fe i. Balmer lines are absent. Along with the earlier results, these features show the star to be a DQZ WD, i.e., having a He-dominated atmosphere, but also containing carbon (Q) and heavier metals (Z). Model-atmosphere fitting to the STIS spectra resulted in a significantly lower ${T}_{{\rm{eff}}}$ of 7740 ± 50 K, and a correspondingly larger radius of 0.01234 ± 0.00032 R_⊙, based on a V magnitude of 10.82 ± 0.03 obtained from the observed SED. This radius, along with a lowered astrometric mass of 0.602 ± 0.015 M_⊙ from G00, removed the earlier discrepancy with the mass–radius relation for CO-core WDs.

Our new dynamical mass for Procyon B allows refinement of a number of its astrophysical parameters and a stringent test of theoretical WD models. We have slightly modified the radius determined by P02, by adjusting for our adopted parallax, obtaining R_B = 0.01232 ± 0.00032 R_⊙. In Figure 9, we compare theoretical predictions with our new parameters for Procyon B. We use theoretical modeling data from the Montreal photometric tables¹⁸ for WDs with pure-helium atmospheres and CO cores. The top panel in Figure 9 shows the location of Procyon B in the theoretical HRD (log L/L_⊙ versus $\mathrm{log}{T}_{{\rm{eff}}}$), along with the model cooling tracks for DB WDs with masses of 0.5, 0.6, and 0.7 M_⊙. The location of Procyon B in the HRD is in excellent agreement with that expected for a WD of our dynamical mass of 0.592 ± 0.006 M_⊙. Also shown in the top panel of Figure 9 are isochrones for ages of 1, 1.25, and 1.5 Gyr, again based on the Montreal tables. By interpolation in the theoretical data, we estimate the cooling age of Procyon B to be 1.37 ± 0.04 Gyr.

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 9, we plot the position of Procyon B in the mass–radius plane. It is compared with a theoretical mass–radius relation for a He-atmosphere CO-core WD with ${T}_{{\rm{eff}}}$ = 7740 K, obtained through interpolation in the Montreal tables. The observed mass and radius are in excellent agreement with the theoretical relation. Also plotted is the Hamada & Salpeter (1961) mass–radius relation for zero-temperature WDs composed of ⁵⁶Fe, which was consistent with the parameters of Procyon B given by Provencal et al. (1997); with our revised parameters, there is no longer agreement with Fe—as first shown by P02.

The surface gravity of Procyon B (in cgs units), based on the mass and radius, is log g = 8.028 ± 0.023. Unfortunately, without Balmer lines, there are no gravity-sensitive features in the HST spectra that would test for consistency with this value. The predicted gravitational redshift is 30.46 ± 0.85 km s⁻¹, but Procyon B possesses no detectable Hα line from which the redshift could be measured using traditional techniques. Onofrio & Wegner (2014) have recently attempted to measure wavelength shifts in the archival HST spectra of Procyon B, using features of Ca ii, Mg ii, and C₂. They appear to have detected the gravitational redshift, but uncertain corrections for pressure shifts are needed.

L13 made a comprehensive analysis of the existing data on Procyon A and B, aiming to determine a consistent picture of the system's evolution. For Procyon B, L13 adopted the P02 effective temperature, but assumed a mass of 0.553 ± 0.015 M_⊙, a value ∼0.04 M_⊙ lower than the dynamical mass determined in the present paper.¹⁹ From these parameters, L13 obtained a cooling age of 1.19 Gyr for Procyon B. Combined with the age of Procyon A, which they determined to be 1.87 Gyr from its position in the HRD, this implied a main-sequence lifetime of only 0.68 Gyr for the progenitor of the WD, corresponding to an initial mass of 2.59 M_⊙. As L13 noted (their Figure 2 and associated text and references), these results placed Procyon B significantly below the initial-to-final-mass relation (IFMR) established from studies of WDs in open clusters. A 2.59 M_⊙ progenitor would be expected to produce a WD with a mass of about 0.69 M_⊙ (cf. Ferrario et al. 2005).

As discussed in Section 7.2, the age of Procyon A may be considerably greater than adopted by L13. This is the case if we use the evolutionary track with large core overshoot, as favored by the GDG14 seismologic analysis. For a Procyon A age of 2.70 Gyr, and our cooling age of 1.37 Gyr for the WD, the main-sequence lifetime of the progenitor of Procyon B was 1.33 Gyr. This corresponds to a ZAMS mass of about 2.2 M_⊙ (if the progenitor had a large core overshoot of β = 1.0 like Procyon A), or a lower initial mass of about 1.9 M_⊙ (if it had a "normal" overshoot of β = 0.2). The Ferrario et al. (2005) IFMR predicts WD masses of 0.65 or 0.625 M_⊙ for such initial masses. Thus there remains a discrepancy, albeit smaller than found by L13, and within the cosmic scatter in the relation (e.g., Figure 2 in L13).

If future observations force an unlikely revision in the current interpretation of the seismic data, Procyon A's age could in principle be as low as ∼1.8 Gyr, resulting in an initial mass possibly as high as about 3 M_⊙ for Procyon B—and a much more severe disagreement with the mean IFMR.

Procyon B presents an unusual case of a WD with the rare DQZ spectral type being a companion of a main-sequence (or slightly evolved) star. A somewhat similar system, HR 637 (GJ 86), was recently studied by Farihi et al. (2013). The K0 V primary in this binary is orbited in a 15.9-day period by a Jovian planet of perhaps 4.4–4.7 M_Jup, detected through RV measurements (Queloz et al. 2000). The K dwarf also has a WD companion in a more distant orbit, with a period estimated at several hundred years. Farihi et al. used HST/STIS to obtain spectra of the resolved WD companion. They found it to be He-rich, with C₂ absorption bands (spectral type DQ6), and having a remarkably similar temperature, mass, and atmospheric carbon content to those of Procyon B. The carbon in cool DQ atmospheres is likely intrinsic to the star. However, the heavier metals (e.g., Ca, Mg, and Fe) seen in Procyon B, but not in HR 637 B, are probably accreted from an external source.

The source of heavy elements accreting onto the photospheres of single DA and DB WDs is usually considered to be a circumstellar debris disk, composed of rocky material (e.g., Jura 2003). Such disks likely form when the WD tidally disrupts terrestrial planets, asteroids, or planetesimals. Although no heavy elements were observed in HR 637 B, Farihi et al. discussed their presence in Procyon B. They examined the heavy-element accretion rates necessary to account for the abundances of Ca and Mg in the Procyon B photosphere, along with the lack of hydrogen, and ruled out accretion from the interstellar medium or a stellar wind from Procyon A. They instead argued for a circumstellar disk around Procyon B as the heavy-element source.

Such a debris disk would be very compact, within a few tenths of a solar radius, and not significantly influenced by the gravity of Procyon A. The object(s) that formed the disk were unlikely to have originated in a protoplanetary disk around the Procyon B precursor, since planet formation would have been confined to within ∼2.3 AU (Holman & Wiegert 1999), and any such planetesimals would likely have been destroyed during the red-giant phases. Instead, Farihi et al. argue that the reservoir could be a much larger disk enclosing the entire binary. They also conclude that it was unlikely that any Jovian-mass objects ever formed around either Procyon A or B, since their formation would have been confined to within the snow limits of each star. This is certainly consistent with the lack of any dynamical evidence for existence of ∼5–10 M_Jup third bodies, as reported in our Section 6.

Another alternative is that the polluting material results from "second-generation" planets, as described by Perets & Kenyon (2013). They suggest that a portion of the wind from an AGB star may be captured by a binary companion (in this case, Procyon A), creating a disk in which planets may form. However, in order to pollute the WD atmosphere, an asteroid or planetesimal born in this disk would then have to be ejected from its orbit around A and into the vicinity of Procyon B.