DISCOVERY OF ROTATIONAL BRAKING IN THE MAGNETIC HELIUM-STRONG STAR SIGMA ORIONIS E

Hesser et al. (1977) observed a primary light minimum of σ Ori E on the night of 1977 January 26/27, as part of their long-term Strömgren uvby monitoring of the star using the number 1 0.4 m telescope at CTIO. Their photometric data were kindly provided to us in electronic form by Professor C. T. Bolton. No error estimates were supplied, so a measurement error Δu = 7 mmag is assumed for all data points—the value quoted by Hesser et al. (1976) as an upper limit on their photometric uncertainties. We do not correct for the apparent season-to-season brightening noted by Hesser et al. (1977),⁶ since this has no effect on the t_min determinations. The u-band light curve is plotted in Figure 1; the accompanying vby data are not shown, since they play no direct role in the period determination (however, see Section 3.4).

Zoom In Zoom Out Reset image size

We observed a primary light minimum on the night of 2004 November 26/27, during Johnson UBVRI monitoring of σ Ori E using the Cassegrain-focus Tek 2048 CCD on the SMARTS 0.9 m telescope (see Oksala & Townsend 2007 for the full data set). For all U-band exposures (in every case, 420 s long), the Tek 1 filter (center 3575 Å, FWHM 600 Å) was used in tandem with the ND3 neutral density filter (7.5 mag attenuation). CCD frames were reduced using standard iraf tasks for bias subtraction and flat-fielding, and cosmic rays were cleaned via Laplacian edge detection (van Dokkum 2001). The phot task of the DAOPHOT package, with an aperture radius of 10 pixels (=396) and a sky annulus radius of 15–20 pixels, was used to perform synthetic aperture photometry. Measurement errors ΔU were estimated as the sum in quadrature of the photometric noise reported by phot, the CCD read noise, and the atmospheric scintillation noise described by the expression on p. 141 of Birney et al. (2006).

As a comparison star, we adopt the nearby (30'') σ Ori D (HD 37468D; B2V; V = 6.62) due to its color and brightness similarity; the U_E − U_D differential light curve is plotted in Figure 1. To make an independent check on the constancy of σ Ori D, we compare its U-band data against σ Ori AB (HD 37468; O9V+B0.5V; V = 3.80). Over the night, we find a mean difference $\overline{U_{\rm D}-U_{\rm AB}} = 3.240\ {\rm mag}$ and a standard deviation σ(U_D − U_AB) = 2.8 mmag; the latter is consistent with the mean measurement error $\overline{\Delta (U_{\rm D}-U_{\rm AB})} = 2.6\ {\rm mmag}$.

We observed a primary light minimum on the night of 2006 January 30/31. The one modification to the 0.9 m instrumental setup was the removal of the ND3 neutral density filter from the light path, thereby shortening exposure times to 2 s (with the idea of obtaining more measurements during the night). In hindsight, this was a counterproductive move. Although the proximity of σ Ori D and σ Ori E mean that shutter corrections are unimportant, the exposures are strongly affected by scintillation noise, especially toward the end of the night when the airmass becomes large. Moreover, σ Ori AB is saturated in all CCD frames and cannot be used as a check star.

The observations were reduced and analyzed as before; the U_D − U_E light curve is plotted in Figure 1.

We observed a primary light minimum on the night of 2008 November 25/26. To obtain more-reasonable (22–65 s) exposure times than in the 2006 season, we reintroduced a neutral density filter (ND2; 5.0 mag). Changes to the default telescope configuration led to the inadvertent substitution of the Tek 2 U-band filter (center 3570 Å, FWHM 660 Å) in place of the Tek 1 filter used previously, but this should have negligible impact on our results. Due to a further oversight, the CCD gain was adjusted from the previous setting of 3.1 e⁻/ADU, to 0.6 e⁻/ADU; this had the unfortunate effect of significantly elevating the readout noise.

The observations were reduced and analyzed as before; the U_E − U_D light curve is plotted in Figure 1. Comparison of σ Ori D against σ Ori AB reveals a mean $\overline{U_{\rm D}-U_{\rm AB}} = 3.236\ {\rm mag}$ and a standard deviation σ(U_D − U_AB) = 12 mmag. The latter value is rather larger than the mean measurement error $\overline{\Delta (U_{\rm D}-U_{\rm AB})} = 9.8\ {\rm mmag}$, but not overly so.

We observed a primary light minimum on the night of 2009 November 28/29. We kept the ND2 filter from the 2008 season, but reverted back to the original Tek 1 U-band filter, and set the CCD gain to 1.5 e⁻/ADU; this led to exposure times between 130 and 250 s. The observations were reduced and analyzed as before; the U_E − U_D light curve is plotted in Figure 1. Comparison of σ Ori D against σ Ori AB reveals a mean $\overline{U_{\rm D}-U_{\rm AB}} = 3.224\ {\rm mag}$ and a standard deviation σ(U_D − U_AB) = 4.2 mmag, the latter consistent with the mean measurement error $\overline{\Delta (U_{\rm D}-U_{\rm AB})} = 4.0\ {\rm mmag}$.

The determination of the time of light minimum (or maximum) in astronomical objects has received significant attention in the literature (see Sterken 2005, and references therein). Historically, the Kwee & van Woerden (1956) algorithm had long been the standard tool, but more recently polynomial fitting has emerged as a powerful and robust approach. Thus, we determine the time of light minimum, t_min, for each of the five light curves plotted in Figure 1 by an adaptive parabola-fitting procedure.

• 1.  
  The time of the dimmest point in the light curve is chosen as the initial t_min.
• 2.  
  Weighted χ² minimization is used to fit a parabola to those points lying within 2 hr of the current t_min (this time interval is chosen to exclude the non-eclipse parts of the light curve).
• 3.  
  A new t_min is chosen by analytically evaluating the minimum of the fitted parabola.
• 4.  
  Steps (2) and (3) are repeated until the value of t_min no longer changes. (Typically, 3–4 of these iterations are required.)

The parabolic fits are drawn over the light curves in the figure, and Table 1 documents the t_min and reduced χ² values associated with each fit. Here and throughout, times are expressed as heliocentric Julian dates (HJD). The tabulated t_min error bounds are calculated from bootstrap Monte Carlo simulations (Press et al. 1992). Specifically, for a given light curve comprising n points, a synthetic light curve is generated by selecting n points at random with replacement. The parabola-fitting procedure is then used to determine t_min for the synthetic curve. This sequence is repeated many times to build up a population of synthetic t_min values, from which we derive a probability distribution function (PDF), F(t_min), reflecting a best estimate of the distribution from which the actual t_min measurement is drawn.

Figure 2 plots the PDFs associated with the five t_min measurements; each is based on a population of 10⁷ synthetic light curves. The error bounds quoted in Table 1 are the confidence limits relative to t_min that each enclose 34.1% of the area under the PDF. Although this choice mirrors the 1σ limits of a Gaussian distribution, the PDFs in Figure 2 underscore that the analysis here does not presume Gaussian errors. Moreover, while the final error in t_min is insensitive to the estimates made for the measurement errors on the individual photometric data, the Monte Carlo simulations used to derive the PDFs do properly allow the actual errors of the noisier data, e.g., in 2008, to produce a larger uncertainty in t_min.

Zoom In Zoom Out Reset image size

We apply the standard O − C diagram technique (e.g., Sterken 2005) to assess the evolution of the rotation period. With the reference epoch defined by Hesser et al. (1976) and the rotation period measured by Hesser et al. (1977), an observed-minus-calculated value

$$\begin{equation} O-C= t_{\rm min}- (2442778\mbox{$.\!\!^{\mathrm d}$}819 + 1\mbox{$.\!\!^{\mathrm d}$}19801 E) \end{equation} \tag{ 1 }$$

is evaluated for each t_min in Table 1; here and in the table, the cycle number E is the integer that minimizes |O − C|.

Figure 3 plots O − C as a function of E; the error bars on each point are taken from the corresponding error bounds on t_min. Also shown in this O − C diagram are fitted linear and quadratic models of the form

$$\begin{equation} (O-C)^{\rm mod}= b_{1} + b_{2} E \end{equation} \tag{ 2 }$$

and

$$\begin{equation} (O-C)^{\rm mod}= b_{1} + b_{2} E+ b_{3} E^{2}, \end{equation} \tag{ 3 }$$

together with the associated fit residuals. The linear model corresponds to a constant rotation period, while the quadratic one represents a period that increases linearly with time. The coefficients {b_j} are determined by a non-χ² maximum likelihood estimation, which seeks to minimize the statistic

$$\begin{equation} Q= \sum - \log F[O-C- (O-C)^{\rm mod}]; \end{equation} \tag{ 4 }$$

here, the summation is over the five light minima, each with its respective O − C, (O − C)^mod, and F(t_min). For a given model, the quantity exp(−Q) is proportional to the likelihood that the measured O − C values could have arisen by chance-fluctuation departures from the model. We use a downhill simplex algorithm (Press et al. 1992) to minimize Q, iterated until the maximum relative deviation between simplex vertices drops below 10⁻⁴.

Zoom In Zoom Out Reset image size

Table 2 summarizes the coefficients {b_j} of the linear and quadratic models. As in Section 3.1, error bounds are determined via Monte Carlo simulations. However, rather than generating synthetic O − C data by bootstrapping, we construct them by perturbing each t_min with random deviates drawn from the appropriate PDF (see Figure 2). From the same simulations we also determine the distribution of the Q statistic, enabling us to associate the linear and quadratic fits in Figure 3 with a likelihood p_null (also specified in Table 2) that the null hypothesis is true: that is, that the deviation of the O − C values from the model arises purely due to chance fluctuations.

Table 2 indicates that the null hypothesis for the linear O − C model is extremely unlikely (p_null = 0.05%), allowing a constant rotation period to be ruled out with a high degree of confidence. However, the converse is true for the quadratic model, which fits the data extremely well. Combining its coefficients with Equations (1) and (3) gives a revised ephemeris for the primary light minimum of σ Ori E as

$$\begin{eqnarray} t_{\rm min}&=& 2442778.8290^{+0.0012}_{-0.0014} + 1.1908229^{+0.0000012}_{-0.0000010} E\nonumber\\ &&+ 1.44^{+0.10}_{-0.11} \times 10^{-9} E^2\ {\rm days}. \end{eqnarray} \tag{ 5 }$$

By taking the derivative with respect to cycle number, the instantaneous period is found as

$$\begin{equation} P = 1.1908229^{+0.0000012}_{-0.0000010} + 2.89^{+0.19}_{-0.22} \times 10^{-9} E\ {\rm days}, \end{equation} \tag{ 6 }$$

which grows linearly at a rate $\dot{P} = 77\ {\rm ms}$ per year. (We discuss the implicit assumption of smooth period growth in Section 4.) At the reference epoch E = 0, this period is rather larger than the P = 1.19081 ± 0.00001 days reported by Hesser et al. (1977); however, these authors' error estimate seems overly optimistic. Reiners et al. (2000) determined a period P = 1.19084 ± 0.00001 days by combining new helium-line data with historical measurements from Pedersen & Thomsen (1977); their value is in good agreement with the mean period $\overline{P} = 1.190833$ days obtained by averaging Equation (6) over the 1976–1998 interval spanned by the helium data.

Before discussing the significance of the measured period increase, we briefly review factors that may have a systematic effect on this result. As demonstrated by Townsend (2008), the timing of light minima is sensitive to the optical depth of the magnetospheric plasma clouds. In principle, the progressive shift in t_min toward later times (seen in the residuals plot of Figure 3) could be explained by the secular accumulation of plasma in the magnetosphere (see, e.g., Townsend & Owocki 2005). However, based on the examples given by Townsend (2008, his Section 2.4), the observed shift of ≈0.02 days between the 2004 and 2009 seasons would require a factor-six increase in the optical depth of the magnetosphere. The minima in Figure 1 clearly do not exhibit this kind of dramatic change. Indeed, although there are some season-to-season variations in the minima depths (on the order of 20 mmag), no monotonic trend is seen; we believe the variations are probably due to changes in the distribution of scattered light from σ Ori AB. Accordingly, we rule out the possibility that the t_min shift is due to plasma accumulation.

Similar reasoning can be addressed to concerns over the change in filters, from Strömgren u in 1977 to Johnson U in the later observations. In passbands where the magnetosphere is more opaque, the time of primary light minimum will tend to occur later. This color–t_min correlation can be clearly seen in the Hesser et al. (1977) observations; in the Strömgren u band (falling blueward of the Balmer jump), the time of primary light minimum is 0.007 days later than in the v band (falling redward of the Balmer jump). Because the U band straddles the Balmer jump, the expected time lag between u and U should be around half of this, i.e., ∼0.003 days. An adjustment of this order to the 1977 t_min point, to correct for the color–t_min correlation, has a negligible effect on our results.

A final possible issue comes from the use of parabolae to measure the primary minimum times. As discussed by Sterken (2005), quadratic fitting is often eschewed in light curve analysis on the grounds that it is unable to adequately model asymmetric light minima. However, if the shape of the light curve does not vary from cycle to cycle, then this bias is irrelevant: it does not matter that t_min occurs slightly before or after the precise time of minimum, as long as the lead or lag remains invariant.

Nevertheless, to explore any bias introduced by the use of quadratic fitting, we have repeated our analysis using cubic fitting to measure the t_min values. Three salient points stand out from this re-analysis: (1) the χ²_red values of the cubic fits are not significantly smaller than those of the quadratic fits (see Table 1); (2) there is no evidence for a systematic lag or lead between the quadratic and cubic minima; (3) there is an obvious difference between the widths of the PDFs, which are a factor of ∼2 broader in the cubic cases than in the quadratic ones. Thus, we conclude that quadratic minimum fitting does not introduce any appreciable bias, and moreover is the more robust approach.