BINARY CENTRAL STARS OF PLANETARY NEBULAE DISCOVERED THROUGH PHOTOMETRIC VARIABILITY. II. MODELING THE CENTRAL STARS OF NGC 6026 AND NGC 6337

The PN NGC 6026 has a lopsided appearance in images taken with short exposure times (see Figure 2). However, deep exposures reveal a symmetric elliptical morphology. Frew (2008) gives a nebular age of 6 × 10³ years and CS temperature of >35,000 K. Stanghellini et al. (1993) give a Zanstra temperature for the CS of 68,000 K and Grewing & Neri (1990) give a temperature of 26,000 K from UV color indices.

Zoom In Zoom Out Reset image size

The original photometric data in Hillwig (2004) were fit using a sine curve resulting in a period with a value half that of the eventually determined orbital period. The new V-band light curve of NGC 6026 folded on this period shows obvious departure from a sine curve as well as significant scatter at minimum and maximum. The solution to these issues is found in the radial-velocity data which show a period of twice that found originally from photometry. The radial-velocity data were used in periodogram, a program which uses the period search technique of Scargle (1982) as modified by Horne & Baliunas (1986). The highest peak is located at P = 0.528 days.

Radial velocities obtained from the spectra as described above were used with the multi-epoch light curve data to determine a more accurate ephemeris. The derived ephemeris for NGC 6026 is

$$\begin{equation} T=2453482.7916(15)+0.528086(4)E, \end{equation} \tag{ 1 }$$

which is in the form T = T₀ + P_orb × E, where T gives the HJD of orbital phase zero (taken to occur at superior conjunction of the CS) and E is the number of elapsed cycles since T₀. Figures 3 and 4 show the radial-velocity data and light curves folded on the ephemeris given above. The fit lines in both figures are from the Wilson–Devinney modeling described below for $q=M_2/M_{\rm{CS}}=1.0$, where M₂ is the mass of the companion and $M_{\rm{CS}}$ is the mass of the CS.

Zoom In Zoom Out Reset image size

Folded on the spectroscopic orbital period, the light curve of NGC 6026 (see Figure 4) shows two photometric maxima per orbit, suggesting that the photometric variation is caused by ellipsoidal variability. The two-cycle-per-orbit light curve and the low amplitude of the variability also indicate that the irradiation effect is not strong. If an irradiation effect is present in a light curve dominated by ellipsoidal variability, it would create unequal depths in the two minima. Minima depths are also affected by the difference in gravity darkening between the inner hemisphere (near the L1 point) and outer hemisphere. In the case of irradiation in an ellipsoidal variability dominated system, the depths of the two minima are then dependent upon both any irradiation effect and gravity darkening in the atmosphere.

We do find deviations in the light curve that do not show periodic behavior. In particular, we see in our phase-folded light curve from phase 0.2–0.4 the 2005 April 17 data are approximately 0.03 mag fainter than overlapping photometry from nights in 2002, 2004, and 2005. The total extent of the data from 2005 April 17 is from about phase 0.9–0.4. Between phases 0.9 and 0.2, the photometry agrees with the rest of the data.

Two other locations in which the light curve appears inconsistent are just after phase 0.6, where data from 2002 April 30 are also approximately 0.03 mag fainter than the remainder of the data. Also, around phase 0.8, data from 2005 April 15 show a similar drop, though in this case lasting for less 10% of the phase. Both before and after, the data correspond to that from the other nights.

There is no immediately recognizable instrumental or atmospheric cause for these drops in brightness. The seeing, transparency, and focus are all stable and the airmass does not drop to unreasonable values (nor would it be expected to cause a fast drop to values that are consistently fainter than the unaffected data for up to 2.5 hr). This suggests that the drop may be intrinsic to the system. In such a case, intervening knots of dust in the surrounding nebulosity are one possible cause. The relatively broad emission lines of N iii and O ii near λ4640 and He ii λ4686 (see Figure 1) are reminiscent of wind emission observed in Of-type stars. However, the gradient in intensity of the nebula at the location of the CS (see Figure 2) makes complete nebular subtraction difficult. Therefore, the He ii λ4686 and Hγ λ4340 lines also include emission due to incomplete subtraction of the nebular spectrum. The N iii and O ii lines near λ4640 also show variable line strength and velocities that do not appear to be associated with the orbital period. The variability also suggests that these lines may originate in a wind, similar to colliding winds observed in other hot star binaries.

The fainter episodes are not just confined to the V-band data, but also appear in the B and R filters as well. Figure 4 shows the phase-folded light curves in all three filters. The B and R filters are more limited in data points because observations in these filters were only obtained during the 2005 observing run.

Zoom In Zoom Out Reset image size

The CS spectrum, apart from emission filling in the Balmer and He ii λ4686 lines, is a good match to HD 47839 an O7 V star. A representative spectrum of HD 47839 is shown in Figure 1 with the CS spectrum. Comparing our spectrum to a tlusty model atmosphere grid (Hubeny & Lanz 1995) produces a good fit at T = 37,500 ± 2500 K and log g = 4.0 ± 0.2, based on the He i λ4471/He ii λ4541 line ratio as well as the wings of the Balmer lines in our spectral range. Our results agree well with the available literature values discussed above.

A Wilson–Devinney code (Wilson 1990; Wilson & Devinney 1971) was used with our photometric and radial-velocity data to find a best-fit solution for the system parameters. For this fitting, we omitted the points described above that fall below the majority of the data. Omitting these data includes the implicit assumption that these data are not representative of the true system brightness, which may or may not be correct. However, these points represent discrete drops in the light curve, so we have used the data which exhibit smooth, continuous light variations. The smoothly variable data allow the model fitting to provide a solution.

There are several observations that can be made to constrain the initial fit parameters. First, the period is quite well determined by the combination of three epochs of photometry and radial-velocity data scattered over roughly the same time period. Also, the fairly short orbital period places some limits on the physical size of both stars. The single-lined radial-velocity curve combined with the orbital period provides a direct relationship between the system mass ratio, $q=M_2/M_{\rm{CS}}$, and the average orbital separation, a. This comes about because the velocity amplitude for the visible CS combined with the orbital period provides a known value for the semimajor axis of the CS orbit, $a_{\rm{CS}}$. Then, we have that $a = a_{\rm{CS}}+a_2=a_{\rm{CS}}(1+q^{-1})$. Likewise, the sum of the masses is related to P and a and thus to P and q. However, since q is the ratio of the masses, then a specific value of q results in unique values for $M_{\rm{CS}}$ and M₂. For the purposes of our modeling, we will assume that any WD or pre-WD must have a mass that falls in the range 0.3 M_☉ ⩽ M ⩽ 1.4 M_☉.

Apart from the above parameters, we also make some specific deductions based upon the light curves. As mentioned previously, the double-peaked curves likely indicate an ellipsoidal effect and because we only see a single-lined spectrum, it is likely that the subdwarf CS nearly fills its Roche lobe. Also, due to the small radius required by the short orbital period, a main-sequence companion would be limited to the lower end of the mass range. The large temperature difference between the two components, combined with the large sizes of the stars relative to their separations, would result in a correspondingly large reflection effect. It would then be the reflection effect that would dominate the light curve.

The light curve also shows that the two maxima occur at different brightnesses. This only happens in ellipsoidal variability for a non-circular orbit. The calculated models show an eccentricity range of 0.010 and 0.025. Furthermore, a close inspection of Figure 4 shows that near phase 0.5 a clear discontinuity in the light curve is present that suggests an eclipse. The possible eclipse also appears to be nearly flat at mid-eclipse suggesting a totally eclipsing system. In fact, the orbital phasing we have used identifies this as an eclipse of the companion, not the CS. There is no clear eclipse of the CS at phase 0.0. Combining these two features suggests that the companion star is small (possible total eclipse, width of the eclipse, and no visible primary eclipse) and that the companion star is very hot (depth of the eclipse and no visible primary eclipse). The model suggests that the most likely companion to the CS is a hot WD or subdwarf. The best-fit model, shown in Table 3, results in a very high surface temperature for the companion. However, the temperature value is directly related to the stellar radius. The two parameters define the total luminosity of the star which directly impacts the amplitude of the variations as well as the depth of the observed eclipse. A small change in temperature requires a corresponding change in stellar diameter. Therefore, the two are strongly coupled in the model fitting. The flat eclipse removes the radius/temperature degeneracy, but only in the presence of a known mass ratio.

Table 3. Best-fit Model Physical Parameters for the CS of NGC 6026

Parameter		Model Values
q	0.63	1.0	1.35
TCS (×103 K)	22	38 ± 3	38
T2 (×103 K)	109	146 ± 15	146
MCS (M☉)	1.40	0.57 ± 0.05	0.30
M2 (M☉)	0.88	0.57 ± 0.05	0.41
RCS (Pole, R☉)	1.41	1.06 ± 0.05	0.78
R2 (R☉)	0.069	0.047 ± 0.005	0.043
$\lg g_{\rm{CS}}$ (cgs)	4.24	4.15 ± 0.02	4.08
$\lg g_2$ (cgs)	6.70	6.85 ± 0.05	6.78
LCS (L☉)	411	2070 ± 80	1121
L2 (L☉)	593	887 ± 40	742
e	0.010	0.025 ± 0.005	0.025
ω (rad)	3.84	3.84 ± 0.9	3.84
i (deg)	80	82 ± 2	84
a (R☉)	2.447	2.87 ± 0.10	3.618

Download table as:  ASCIITypeset image

Using the mass limits for the subdwarf/WD companion described above, with the relationship between q and the masses, provides a mass ratio range of 0.63 ⩽ q ⩽ 1.35 and a corresponding range in CS mass of $1.40\,M_\odot \ \ge M_{\rm{CS}}\ge 0.30$ M_☉. Varying the mass ratio affects the stellar separation, a, the stellar radii, and to a lesser extent the surface gravities. For a given mass ratio, increasing the CS mass in turn directly affects the stellar separation, the required size of the CS (to achieve the ellipsoidal effect), and thus the required size (and temperature) of the companion. However, once again, the suspected eclipse provides some limits to these parameters. The relative stellar radii are constrained by three items: the width of the eclipse, the subdwarf's Roche lobe filling factor required for the ellipsoidal effect, and the relative amplitude of the eclipse and ellipsoidal effect. Putting all of this together means that for a given mass ratio, we derived reasonably constrained physical parameters for the system.

In Table 3, we show the results of our modeling for a mass ratio, q = 1. The errors given in the table represent the spread of values over which each variable still produces a good fit. We also show in the table the range for each value resulting from varying the mass ratio between the two limits discussed above. These ranges represent the possible values based on physical limits in the system. Mass ratios approaching the lower limit correspond to the largest component masses because the radial-velocity amplitude of the CS is set by the data. For the low q model the CS temperature is forced to values too low to be associated with our observed spectrum. To maintain a reasonable temperature for the CS, the mass ratio range can effectively be limited beyond that shown in Table 3, with an approximate lower limit of 0.8 M_☉ ($M_{\rm{CS}}\approx 0.83$ and M₂ ≈ 0.67). Furthermore, we may also expect that if the companion is more compact than the CS, then it is evolutionally more advanced than the CS. For a more evolutionary advanced state, it must have had a more massive progenitor and would likely have a more massive remnant. In this case, the mass ratio would be q ≳ 1.0. The modeled eccentricity, temperature, and inclination remain the same for mass ratios in the range 1.0 ≲ q ≲ 1.35. For values below q = 1.0, each of these physical parameters must be changed to result in a good fit to the data. The best-fit values shown in Table 3 correspond to the broadest mass ratio limits, q = 0.63 and q = 1.35. Also shown are model values for the likely lower limit of q = 1.00. We expect that the system parameters lie between the models for q = 1.00 and q = 1.35. Approximate uncertainties for individual parameters are given in the q = 1.00 model.

The range of values for the stellar diameter of the companion between 0.043 and 0.069 R_☉ places it at more than twice the expected radius of an object supported by electron degeneracy pressure (see, e.g., Schmidt 1996). Both the larger than expected radius and high temperature of the companion may be due to accretion heating via a CE phase or an AGB/RGB wind. The accreted material may have created an expanded hot atmosphere. It does not seem likely that such an expanded atmosphere due to CE accretion would still be present given the approximate 6000 year age of the PN (Frew 2008) and the cooling rates for WDs in novae which occur on the scale of hundreds of years (Prialnik 1986). It is possible, given the non-cyclical variations seen in the light curve described above and the evidence in the spectrum for winds, that more recent wind accretion may be responsible for the increased size and temperature. It is unlikely that the faint, hot companion would be younger than the larger CS, given its much smaller size. It is possible that the companion has not yet contracted to the WD cooling track, which would require that the masses of the two progenitor stars were very nearly identical (the q = 1.0 model), which is statistically very unlikely.

We have used passband specific linear limb darkening coefficients from Van Hamme (1993) for the models. For the compact companion, the limb darkening coefficient used was 0.069 in all three filters. The CS coefficients used were 0.256 (B filter), 0.199 (V filter), and 0.095 (R filter). The final result is sensitive to the choice of limb darkening coefficient only for the CS and only for changes of the order of 20% or more in the coefficient value. Our values should be well within this uncertainty. In addition, given the high temperatures for both stars, we used an albedo of 1.0 for both components.

Comparing our temperature and log g values to the model atmosphere analysis of Napiwotzki (1999) in his Figure 2, we find that the CS of NGC 6026 falls very close to the evolutionary track for a 0.546 M_☉ post-AGB star and a 0.414 M_☉ post-RGB star. Our range in masses covers both of these values, so that we are unable to distinguish between a post-AGB versus post-RGB evolution for the CS based on the model results. However, the nebula is very highly ionized (de Vaucouleurs 1955), lacking [O ii] line emission and with strong He ii emission. The high ionization is inconsistent with the previously observed 35,000 K CS, but is more in line with the >109,000 K companion discussed above. It then appears that much of the ionization of the nebula is due to the compact, hot companion rather than the visible CS.

The modeled systemic velocity of the binary is γ = −93.4 ± 2.1 km s⁻¹ which is consistent with the heliocentric radial-velocity for the PN in the literature of γ = −101.5 ± 4.7 km s⁻¹(Durand et al. 1998).

NGC 6337 appears to be a nearly perfectly circular PN in most images (see Figure 5) and has been classified in most lists as round (or elliptical in the case where round nebulae are not classified separately). However, Corradi et al. (2000) and García-Díaz et al. (2009) have shown that the round nebulosity seen is the narrow waist of a nearly pole-on bipolar PN. The papers present deep images and spectra showing the edges of the polar lobes extending slightly outward from the bright ring in the northwest and southeast directions with detected polar caps. The nearly circular appearance of the bright ring suggests that the major axis of the nebula is tilted by ≲15° to our line of sight. The distinct appearance of the lobes on opposite sides of the ring, however, suggests that the tilt is likely closer to 15° than 0°. The true inclination though is difficult to extract without some knowledge of the axial ratio of the nebulosity, which can vary drastically among bipolar PNe.

Zoom In Zoom Out Reset image size

It is also interesting to note that Corradi et al. (2000) and García-Díaz et al. (2009) also detect point-symmetric structures suggesting precessing outflows in their deep images of NGC 6337. The structures are reminiscent of similar features seen in the lobes of other bipolar PNe (e.g., Sahai 2004). Also seen in García-Díaz et al. (2009) and in our Figure 5 as the red filaments are features reminiscent of FLIERS (Balick et al. 1993).

In the case of NGC 6337, there are several published values for physical parameters of the nebula and CS. Regarding the temperature of the CS, the results of Gorny et al. (1997) and Stanghellini et al. (1993), 41,700 K and 44,870 K, respectively, are in good agreement and are both H i Zanstra temperatures. Amnuel et al. (1985), however, report a CS temperature of 105,000 K using nebular diameter and 5 GHz radio flux. The CS temperature is discussed below in the context of model fits to our data.

Gorny et al. (1997) present, in addition to CS temperature, estimates of the CS mass, evolutionary time of the nebula, and distance to the PN. Their distance, 1.6 kpc, agrees well with the published values of Acker et al. (1992, 1.1–1.7 kpc) and Stanghellini et al. (1993, 1.32 kpc). However, Frew (2008) gives a distance of 860 ± 200 pc. The CS mass of 0.563 M_☉ given by Gorny et al. (1997) is consistent with their relatively low temperature. Evolutionary models then provide their reported evolutionary time of 2.42 × 10⁴ years. Frew (2008) reports an age of 1.2 × 10⁴ years derived from his surface brightness/radius relation. The major difference in nebular ages likely results from the author's different distance values.

As with NGC 6026, the CS of NGC 6337 was reported to be a likely close binary by Hillwig (2004). The preliminary orbital period was 0.1735(3) days. The ring-like appearance of NGC 6337, with a short period binary CS, is reminiscent of the PN Sp 1 and its close-binary CS, which has an orbital period of 2.9 days (Bond & Livio 1990).

The orbital period given by Hillwig (2004) was determined using a sine curve fit to the photometry. Our additional photometry shows that the light curve of the CS continues to be well fit by a sine curve and that there is no apparent ellipsoidal effect as with the CS of NGC 6026. The periodogram using all of the available photometry shows strong aliases relating to the time gap between observation epochs. The period given below provides a dispersion in the data which is nearly half of the next best period and has been taken to be the correct orbital period. The derived photometric ephemeris (in HJD) for minimum light in the CS of NGC 6337 is

$$\begin{equation} T=2452395.5368(1)+0.1734742(5)E. \end{equation} \tag{ 2 }$$

The near sine-wave appearance of the phase-folded light curve along with the short orbital period suggests that the cause of variability is an irradiation effect. Effective modeling of such an effect is nearly impossible with only one filter due to the high correlation between many of the important physical parameters. The multi-color photometry allows a more accurate determination of the system's physical parameters, though lack of a radial-velocity curve still limits the precision of the fitting. A more thorough description of these issues is given in Paper I.

We have once again used the Wilson–Devinney code to model the available data. Immediately apparent in Figure 6 are the very similar amplitudes in all three filters. The CS here is surrounded by a faint, very uniform nebulosity (inside the bright ring), so light contribution from the nebula should be negligible. The small difference in amplitude between the three filters then provides a very useful discriminating tool for the model fitting. Initial models were calculated using $M_{\rm{CS}}=0.6$ M_☉ and $T_{\rm{CS}}=45{,}000$ K, the approximate values from Gorny et al. (1997) and Stanghellini et al. (1993).

Zoom In Zoom Out Reset image size

The short orbital period also provides some useful constraints on model fitting. The required small separation means that a secondary star with M₂ ≳ 0.35 M_☉ would have a main-sequence radius larger than its corresponding Roche lobe, though most close-binary CSPNe with irradiation effects appear to have secondary stars with radii greater than their main-sequence values (Paper I) which would require a lower value for the secondary mass upper limit. The well-behaved variability we see during three epochs separated by about three years suggests that the secondary star is not in contact with its Roche lobe. Therefore, we have an upper limit for M₂ in our models. We use this upper limit as a starting point for the secondary mass. Limb darkening and gravity darkening coefficients were determined as for NGC 6026. We use an albedo of 1.0 for the hot CS and 0.5 for the cool companion. One potential problem in models of irradiated binaries is the lack of azimuthal energy transport. The amount of energy transport is not known so models are not currently able to account for it. If such transport occurs, then the amplitudes derived from modeling would be too large as the "dark" side of the companion star would be brighter due to the energy transported from the "day" side.

From the calculated models, we find that our beginning value of $T_{\rm{CS}}$ makes it very difficult to accurately model the observed light curves. In order to arrive at B, V, and R amplitudes that are approximately equal requires that the CS be a subdwarf with a radius equal to a significant fraction of its Roche lobe size. This is not expected from the results of Gorny et al. (1997) which include a long evolutionary time, after which the CS should have contracted to a radius closer to its WD value. We also found it necessary to increase the companion effective temperature to T₂ = 5500 K, much higher than that expected for the 0.35 M_☉ secondary star. However, it is commonly found in irradiated systems that the cool companion has a temperature much higher than that expected for its mass and from what we know of irradiated binary CSPNe, 5500 K is not unreasonable (Paper I).

Figure 6 shows our best-fit curve (dashed line, model 1) for our initial $T_{\rm{CS}}=45{,}000$ K. Also shown is our best-fit curve for an alternate value of $T_{\rm{CS}}$ (solid line, model 2) as discussed below. The parameters for model 1 are given in Table 4.

Table 4. Model Physical Parameters for the CS of NGC 6337

Parameter	Model 1	Model 2
$T_{\rm{CS}}$ (×103 K)	45	105
T2 (×103 K)	5.5	2.3
$M_{\rm{CS}}$ (M☉)	0.60	0.60
M2 (M☉)	0.35	0.20
$R_{\rm{CS}}$ (R☉)	0.09	0.016
R2 (R☉)	0.42	0.34
$\lg g_{\rm{CS}}$ (cgs)	6.30	7.78
$\lg g_2$ (cgs)	4.74	4.67
i (deg)	28	9
a (R☉)	1.286	1.214

Download table as:  ASCIITypeset image

We expanded our modeling to higher $T_{\rm{CS}}$ values in an attempt to obtain a more robust fit to the data. We find that for $T_{\rm{CS}}\gtrsim 90{,}000$ K, we are able to model the variability amplitudes in all three filters effectively over a fairly wide range of initial parameters. For our model 2, then we set $T_{\rm{CS}}=105{,}000$ K, the value from Amnuel et al. (1985). Using the upper limit for M₂ above with T₂ = 3500 K (a reasonable value for an MS star of this mass, see, e.g., Drilling & Landolt 2000, and references therein), we found a strong correlation between the radius of the CS and the system inclination. For i ≳ 10°, the CS radius is greater than that expected for a WD star, as in model 1. For these models, the CS dominates the luminosity in the observed filters. For i ≲ 10°, we find a CS radius of the order of that expected for a WD (roughly 0.01–0.02 R_☉). We also find that the luminosity of the secondary becomes important (due to the much smaller radius of the CS) to the relative amplitudes of the light curve in the three filters.

Rough limits on the orbital inclination, i, can be placed by the models. Variability in close-binary stars due to irradiation generally produces sine-wave like light curves. However, for increasing inclinations, the light curve maximum becomes narrower and the minimum correspondingly wider. Values for i ≳ 40° result in minima that are obviously too broad to fit the data well. Also, good fits for i ≲ 6° are unphysical, even if the CS has fully reached the WD cooling track, as they require $R_{\rm{CS}}$ values significantly smaller than expected for a star supported by electron degeneracy pressure. The inclination of the nebula from Corradi et al. (2000) of ≲15° falls within our modeled range. The close alignment of the binary and nebular axes itself does not unambiguously support a direct physical connection between the CS and the nebular shape, though the probability of a chance alignment is very low. In addition, it is likely that the rotation axis of the CS is aligned with the orbital axis. So while the alignment of binary axis and nebular axis is evidence for a direct physical connection, it does not allow us to differentiate between shaping models.

As mentioned above, the secondary star becomes important for inclinations below about 10°. In all cases, the secondary star very nearly fills its Roche lobe and the values of T₂ decrease from 3500 K at 10°–11° to about 1500 K at 6°–7°. Figure 6 shows the best-fit model at 9° (model 2, solid line) and the corresponding parameters are given in Table 4. Below about 9°, the secondary temperatures are probably unrealistic, again given the higher than normal temperatures seen in irradiated companions (Paper I).

One of the clearly distinguishing characteristics between models which fit the data lies in the radial-velocity curves of the components. The range in inclinations produces models with radial-velocity amplitudes which vary drastically. Fortunately in this case, it is likely that a spectroscopic study could result in a double-lined radial-velocity curve and we are pursuing these observations. The double-lined spectrum is due to the observed nature of the irradiation effect. The CS, which dominates the system luminosity in the optical range, may show absorption lines, while the irradiated hemisphere of the secondary star will show strong emission lines. A spectroscopic program should also be able to distinguish between the high and low temperature models by using non-local thermodynamic equilibrium codes to model the absorption lines of the CS. The limitations on the system's inclination provided by radial velocities would also provide stronger constraints on the likelihood of a physical connection between the inclination of the binary CS and the inclination of the nebula.

The modeling results for a 0.6 M_☉ CS show that the photometric period is likely to be the same as the orbital period and due to an irradiation effect from the CS on a cool companion of mass M₂ ≲ 0.35 M_☉. The system inclination appears to be ≲40° but is more likely 6° ≲ i ≲ 20°. The broad ranges of reported CS temperatures in the literature are all supported by the models, though the most robust fitting occurs for $T_{\rm{CS}}\gtrsim 90{,}000$ K. The lack of radial-velocity data means that strong limits cannot be placed on $M_{\rm{CS}}$ or a, the orbital separation. Although the numbers given above apply for one model, given a specific CS mass, they are representative of the range of models for physically appropriate CS masses and provide a general physical description of the system.