Correlated X-Ray and Optical Variability in the O-type Supergiant ζ Puppis

Broadband light curves were constructed for each ObsID, then combined into a single light curve for the entire Chandra campaign consisting of count rates per 4 ks bin. Counts from HEG and MEG were combined where their spectra overlap in wavelength (3.1–10.3 Å). Then the MEG counts in the range 10.3–20.7 Å were added to this total, providing a total wavelength coverage of 3.1–20.7 Å. The break at 10.3 Å was chosen because some of our ObsIDs did not have HEG coverage beyond this position. This break in wavelength coverage also provides compatible ranges to XMM-Newton hard and XMM-Newton medium bands (note that the definitions of "hard" and "medium" differ between XMM-Newton and Chandra; see Table 4). In this paper we will refer to the Chandra definitions of hard, medium, and soft, and an additional band that we call "hybrid-hard" that refers to all of the Chandra medium band and part of the Chandra hard band, in order to directly compare to the XMM-Newton hard band. All times in the light curve were corrected to the barycentric times. The bins are contiguous and exclusive, which means there is always a short time bin at the end of each observation. We binned to 4 ks intervals uniformly from the beginning of each observation. Short end bins (≤50% exposure per bin) were excluded from the timing analysis due to their larger variance, which could degrade results. A total of six short time bins were excluded, containing less than 1.5% of the total exposure.

As in Nazé et al. (2018), we first estimated the amount of variability in the light curves using the indices VI and F_var (Edelson et al. 2002, see Appendix A). VI is a variability index, expressed as $(\max -\min )/(\max +\min )$, with max and min values taken without any exclusion and thus making them prone to noise fluctuations; it thus compares the amplitude of the count rate variation (i.e., half the peak-to-peak variation amplitude) relative to the mean. The more useful fractional variability amplitude F_var provides an idea of the amplitude of intrinsic changes, relative to the mean, as it is defined as

$$\begin{eqnarray*}&&{F}_{\mathrm{var}}=\sqrt{{S}^{2}-{\sigma }_{\mathrm{err}}^{2}}/{X}_{{\rm{m}}},\end{eqnarray*}$$

with the mean X_m = ∑X_i/N, i = bin index, the dispersion ${S}^{2}=\sum {({X}_{{\rm{i}}}-{X}_{{\rm{m}}})}^{2}/(N-1)$ and the mean Poisson error ${\sigma }_{\mathrm{err}}^{2}=\sum {\sigma }_{\mathrm{err},{\rm{i}}}^{2}/N$. This eliminates the variance that is due to Poisson noise (Edelson et al. 2002). A comparison of the variability in Chandra and XMM-Newton data using these indices is presented in Table 2. Note that both Chandra and XMM-Newton values agree with each other and, as found in XMM-Newton data (Nazé et al. 2018), the 10.3–20.7 Å band appears more variable than the 3.1–10. Å band in Chandra data (see Table 4 and Section 4).

The corrected Chandra light curves were analyzed using a modified Fourier algorithm adapted to data sets with uneven sampling (Heck et al. 1985; Gosset et al. 2001; Zechmeister & Kürster 2009).¹⁶ Figure 2 shows the obtained periodogram with a red horizontal line representing the 1% significance level. This level was estimated by two methods, which agree with each other: (1) by performing 2000 Monte-Carlo simulations drawing at random the individual count rates from a Gaussian distribution with the same mean for all simulations and a standard deviation equal to the dispersion of the observed count rates, and (2) by shuffling the data. The maximum amplitudes reached by the periodograms of the simulated data were recorded and the significance level was fixed at the amplitude above which only 1% of the periodogram maxima lie. This 1% significance level thus provides the amplitude that will be detected in a data set similar as ours by chance; some may consider it as an upper limit because the simulations include the effects of real variations (two periodicities along with all their many harmonics and aliases), which will enhance the dispersion of the observed count rates beyond that of the Poisson noise level.

Zoom In Zoom Out Reset image size

Because the sampling of the Chandra observation is very irregular, many aliases are present in the periodogram (see the spectral window in Figure 2). In particular, the data set is composed of short snapshots spread over a year, so any single peak will appear as a broad compound of narrow subpeaks (see zoomed images at the bottom of Figure 2). In the periodogram, two groups have the highest significant frequencies in units of cycles per day (c/d): near 0.17 c/d⁻¹ and near 0.56 c/d⁻¹. Since the second one has the largest amplitude, we first focus on it. A sum of 10 Gaussians was fit to the periodogram over the region of the strongest peaks (lower-right panel in Figure 2). Of the 10 peaks fit, the peak with centroid 0.561147c/d, indicating a period of 1.7820 day, was chosen as the closest to the BRITE period and therefore the most likely true period.

We estimated the error in the period to be of the order of 1/10T = 1/(10 × 411)d = 0.00024c/d, with T being the total time interval of all the Chandra observations in days (as in Ramiaramanantsoa et al. 2020, in preparation). The method yielded a period error of 0.0008 day. We checked this error value using the Gaussian fits described above, where the uncertainty in the centroid value, the sigma value, and the HWHM are consistent with 10% of the peak width and with our adopted error value.

The highest peak at P = 1.7727 ± 0.0008 day within this group of subpeaks is not within the errors of the X-ray period of 1.7820 day and is probably due to the irregular sampling; occasionally in such cases, a subpeak or alias (rather than the highest peak) corresponds to the actual signal. The presence of nearly the same periodicity in both optical and X-ray data sets clearly points toward a common origin for both variability phenomena.

Given that the optical BRITE value is consistent from several years of data (see notably Ramiaramanantsoa et al. 2018), we adopt the BRITE ephemeris (period P = 1.7806 day and T₀ is HJD = 2,458,425.800) for any calculation of light-curve phases related to this period. Using this ephemeris places the optical maximum at phase 0.45, The X-ray light curves are always shown with that ephemeris, in order to compare the light curves directly. Figure 3 shows both X-ray and optical light curves folded with this ephemeris. For the binned X-ray light curve, we calculated averages of data points in 20 phase bins over the 1.7806 day period. Each phase bin is ϕ = 0.05. The BRITE light curve has been binned in the same manner. The peak-to-valley amplitude of this binned light curve is about 6%, which is remarkably large compared with the optical variations (peak-to-valley amplitude of about 1%).

Zoom In Zoom Out Reset image size

The light curves in Figure 3 are notable for several features. The minimum of the X-ray light curve corresponds to the maximum of the optical light curve, approximately. The BRITE light curve has a clear maximum at ϕ = 0.45 as well as a secondary maximum at about ϕ = 0.8. The Chandra light curve has a maximum at  ∼ ϕ = 0.9 and a secondary maximum at about  ∼ ϕ = 0.55.  The phase of the primary X-ray maximum was determined by cross-correlation between the BRITE light curve and the Chandra light curve, yielding an offset of 0.45 in phase. The secondary X-ray maximum is based on the "point" at 0.55 in phase, but this point represents the mean of 18 time bins from 12 observations that fell in this phase bin after folding on the 1.78 day period. Each point in Figure 3 was calculated in the same way, but the number of time bins that fell in each phase bin varies from 5 to 21 due to the uneven coverage of the observations. The data point for the secondary X-ray maximum in Figure 3 is calculated to be somewhat more than 2σ from the mean of the residuals after fitting the light curve with a 2° polynomial. Although this is a marginal detection, it nevertheless should be explored and is discussed in Section 6.5. A nearby point at about ϕ = 0.65 has approximately equal count rate as the peak at ϕ = 0.55, but is less than 1.5σ from the mean of the residuals.

Once we applied the same period and ephemeris to the BRITE data and Chandra data reported here, the peaks of these two light curves did not coincide. The maximum X-ray emission lags by Δ(ϕ) ∼ 0.45 behind the maximum optical emission, a value determined by cross-correlation of the two light curves. Such a large time lag would place constraints on the connections between the optical and X-ray emission regions, as will be discussed in more detail in Section 6. Alternatively, if we take both the primary and secondary maxima in the optical light curve and relate them to the secondary and primary X-ray maxima, we have two optical events with X-ray lag time of  ∼ ϕ = 0.1 each.

When taken as a whole, this Chandra data set for ζ Pup provides coverage of all phases of a 1.78 period. However, not all phases are covered equally well in all observation segments. This can lead to some limitations. For example, an important question is whether the period detections are stable over time. At first glance, since the Chandra data have been mostly obtained in two observing windows, summer 2018 and summer 2019, they could allow for such a check. Unfortunately, as shown by the observation coverage (Figure 4), the maximum of the 1.7806 day modulation was not sampled during the second observing window. However, it is important to note that the binned light curves of the two observing windows, when they do overlap in phase, appear compatible, favoring the hypothesis of the periodogram stability.

Zoom In Zoom Out Reset image size

Figure 5 compares the XMM-Newton calibration observation taken in 2019 April (blue circles) to the Chandra data. The XMM-Newton minimum agrees approximately with the Chandra minimum. The light-curve evolution appears different, with a steeper declining slope for the XMM-Newton data. However, one needs to keep in mind that while the string of XMM-Newton data points in the figure represent a single exposure, the Chandra big dots represent data points at these phases from a combination of the whole campaign. One would not expect exact agreement between XMM-Newton and Chandra because XMM-Newton data are sensitive to softer X-rays than Chandra data, and because of the difference in sensitivity between the XMM-Newton and Chandra instruments. Looking at the small dots of Figure 5, which represent the individual Chandra 4 ks bins, we see that the XMM-Newton binned light curve appears well within the scattered Chandra points.

Zoom In Zoom Out Reset image size

To further inquire about the presence of additional coherent signals or of stochastic variability beyond the 1.7806 day signal, we once again calculated averages of data points, but this time in 10 phase bins rather than 20; the choice of 10 phase bins per cycle is a compromise between having enough signal in each bin while still allowing us to examine the shape of the phased light curve in detail. Linear interpolations of this binned Chandra light curve were used to remove from the Chandra data the variations associated with the 1.7806 day period. Because the light curve does not appear to be a perfect sine wave (see Figure 6), there must be harmonic content in addition to the fundamental. While an improved data cleaning for this period would take into account the contributions of these additional harmonic components, our periodogram (Figure 2) does not provide enough information on the harmonic content to allow for this. Hence the choice of the subtraction of the mean light curve.

Zoom In Zoom Out Reset image size

We then performed a period search on the resulting cleaned light curve (Figure 6). The overall peak amplitudes appear largely reduced in this periodogram, indicating that the 1.7806 day signal dominates the X-ray variability. However, one peak is still significant, at a frequency of 0.16860 ± 0.0002 day⁻¹, corresponding to a period of 5.9312 ± 0.009 day. The next largest peak, only slightly lower significance, lies at 0.1978 ± 0003 day⁻¹, corresponding to a period of 5.056 ± 0.008 day. It is actually difficult to choose which one of those two peaks corresponds to the "real" signal. In fact, different processing choices may lead to one or the other peak leading in the cleaned periodograms. This indicates that, most probably, a period of about 5 days or about 6 days is present.

Focusing on the formally significant period only, the corrected light curve, cleaned for the 1.7806 day signal, was folded with a 5.9312 day period and averages in 10 phase bins were then calculated. The original light curve was then cleaned by this average 5.9312 day signal in the same manner as done before for the 1.7806 day case. The periodogram of this newly cleaned light curve is shown in Figure 7. It is immediately obvious that the 1.7806 day signal remains significant. It is thus important to note that both signals are distinct, i.e., even though they may be close to a harmonic ratio of 3, they appear separate as cleaning by one does not remove the other. As a last trial, we cleaned the original light curve by both averaged curves and performed a period search on the result (Figure 7). This time, the periodogram amplitudes are very much reduced and no significant signal is detected.

Zoom In Zoom Out Reset image size

As a final exercise, we calculated the fractional variability amplitudes F_var of the light curves at different stages of cleaning (Section 3). F_var was 0.031 ± 0.004 for the original light curve; it decreased to 0.020 ± 0.005 after cleaning for the 1.7806 day signal, or 0.024 ± 0.004 after cleaning for the 5.9312 days signal, and finally 0.011 ± 0.007 after cleaning by both signals. The latter value implies that zero (i.e., no variability beyond Poisson noise) is only at 1.2σ. Hence the observed variability of the Chandra light curve can be explained a posteriori by two periodic signals (1.7820 day and 5 or 6 days), with the presence of additional variations, stochastic or coherent, to be confirmed (or excluded) with higher-quality data in the future.

We also plotted in Figure 8 the Chandra light curve folded on a 5.06 day period, which was previously proposed by Moffat & Michaud (1981). This light curve is possibly consistent with the presence of an ∼5 day period, although there are several significant outliers. We discuss this potential period in Section 6.

Zoom In Zoom Out Reset image size

To aid in visualization of the time variability of the emission line fluxes, Figure 10 shows a period-folded, "trailed" spectrum where the vertical axis represents the phase of a 1.7806 day period, using the ephemeris in Table 3. All 9 ks, time-sliced, and calibrated spectra were used in constructing the dynamic image. These images and residuals allow one to view the variations across the spectrum, including emission lines and continuum. Details are in the caption of Figure 10.

Zoom In Zoom Out Reset image size

The Chandra data have modest resolution and significant noise. Line moments of order zero to 2 were thus estimated for strong isolated lines (Si xiv 6.182 Å, Mg xii 8.421 Å, Ne x 12.132 Å, Fe complex near 11.5 Å, and Fe complex near 15.01 Å).¹⁷ This was done in both real and simulated 9 ks slices.

We evaluated the significance of the results of moments by creating a set of "simulated" data using the same time filter and calibration files as the real data, with a constant model corresponding to the mean of the full data set and Poisson noise applied randomly. We then examined these data with the same techniques as the real data. Because they will have the same gross statistical properties and window function, analysis of the simulated data in parallel with the actual data gives us another tool to evaluate the reality of any variability we find. χ² tests against constancy reveal the moments to be in general significantly (i.e., SL < 1%) variable in real spectra, but not in the simulated data set, except for line width. Therefore, stellar variability seems to be present at least in flux and centroid values. When phased with the 1.7806 day period, no coherent flux variations are obvious, but this is probably due to the large noise when considering a single line.

To further assess the variations, we derived the cumulative distribution functions of the moments, both for real and simulated data, and tested their difference using a Kolmogorov–Smirnov test. No significant difference between simulated and real data was found with this test. We also determined the Pearson correlation coefficients between moments of the same lines, again for both simulated and real data. These coefficients never reach values beyond 50% and there is no systematic difference between simulated and real data (e.g., systematically larger coefficients for real data), hence we rule out the presence of significant correlations between the different moments of a line (e.g., flux-centroid, flux-width, and centroid-width correlations).

Finally, the same period-search algorithms were applied to the moments of order zero through 2, for both real and simulated data. No significant difference was found, i.e., peaks with similar amplitudes are found in both data sets. In this case, moment analyses appear too insensitive to allow detection of periodicities in such noisy data.

In the full spectral range, a period of 1.7820 days is clearly identified in the X-ray data and is almost certainly closely connected to the previously measured optical modulation with essentially the same period. The peak-to-valley amplitude in our observations for the P = 1.7820 day signal is about 6%, which is remarkably large compared with the optical variations of only 1%. The 1.7820 day period we found in the Chandra full-band data is not peculiar to some small region of the X-ray spectrum, but is relatively general through the whole range of wavelengths (see Figure 10). It is present in the flux from the H-like lines, He-like lines, and a collection of Fe lines. It even appears to be present in the portions of relatively line-free continuum. The light-curve characteristics, however, appear slightly different, both in shape and amplitude, for different bands or lines. A closer investigation of the emission lines (though strongly limited by the small number of photon counts in each line) did not reveal clear periodic changes in line properties other than flux. Moment analysis of a number of emission lines individually does not indicate any correlation between line width and centroid velocity.

Considering only the maxima in the optical (ϕ = 0.45) and X-ray (ϕ = 0.9) light curves, the lag of ϕ ∼ = 0.45 of a cycle in phase might provide an important clue in untangling the connection between the two bands. Assuming the maximum X-ray flux lags the maximum optical peak is a possible interpretation of the relation between the optical and X-ray light curves due to the amplitudes of the maxima. However, it is also possible that the X-ray maximum is actually physically related to the secondary optical maximum because the two phenomena occur at almost equivalent phases (secondary optical: ϕ = 0.8;  primary X-ray: ϕ = 0.9) and the lag time would be  ∼ ϕ = 0.1. Also, the secondary maximum in the Chandra light curve may be related to the BRITE maximum (primary optical: ϕ = 0.45; secondary X-ray: ϕ = 0.55), again because they are at nearly the same phase. Together, the lag time for this set of features would each be  ∼ ϕ = 0.1.

In addition to the main 1.7820 day X-ray period, there is an indication of some periodicity close to 5 days or to 6 days. This periodicity is more difficult to pin down for two reasons: (1) the breadth and complexity of this peak in the periodogram (see Figure 2, lower-left panel), and (2) the lack of any comparable periodicity in the contemporaneous BRITE data. Though simply choosing the highest peak in this region of the periodogram would favor a period near 6 days, the nearly equally prominent peak near 5 days (see Figure 8) is particularly interesting because of its proximity to two previously measured periodicities for this star. Moffat & Michaud (1981) claimed a 5.1 day periodicity in the near-central reversal in the Hα emission line, and explained it by excess magnetically confined plasma. Howarth et al. (1995) found a 5.2 ± 0.7 day period in two solid uninterrupted weeks of IUE UV spectra, described as being of unknown origin. The light curve of Chandra data folded on a 5.056 day period (Figure 8) suggests that this period is possible, although not as convincingly as the 1.78 day period.

Finally, a residual stochastic component to the X-ray variability cannot be ruled out. After removing the two periodicities described above from the signal, the scatter in the residual broadband X-ray light curve is ∼6%, compared with the estimated Poisson noise of 4%. Because there should not be any significant additional instrumental sources of noise, this leaves a net rms scatter of $\sqrt{({6}^{2}-{4}^{2})}=\sim$ 4.5%, which, if real, would be intrinsic to the wind. However, as shown above from the F_var determination, the presence of stochastic variability is only a 1.2σ detection, hence the need for confirmation. If confirmed, this might be understood in terms of stochastic shocks throughout the wind leading to random fluctuations at this level in the X-ray flux. The X-ray modulations are 0.03 for the initial data set, 0.02 after cleaning for one period, and 0.01 after cleaning by both periods, i.e., the 1.78 day and 6 day periods may have similar strengths as the remaining variability. A similar ratio of periodic versus stochastic variability was also found in the optical (Ramiaramanantsoa et al. 2018). Nevertheless, because of the Poisson noise, additional data are clearly needed to clarify the presence of this stochastic variability.

Assuming for the moment that this stochastic component exists and will be confirmed in the future, one is thus drawn to the idea that we are seeing similar effects in X-rays as in optical data regarding CIRs and wind clumps. If the X-ray emission showed no intrinsic stochastic emission, it would indicate a scenario of Poisson saturation due to myriad clumps. Nazé et al. (2013) found that the ζ Pup wind should contain at least 10⁵ X-ray-emitting shocks, leaving a relative Poisson fluctuation of $\sim 1/\sqrt{{10}^{5}}=$ 0.3%, well below detection limits of current X-ray telescopes. However, the modeling for this estimate neglects the more likely scenario of a turbulent power law with progressively fewer clumps or shocks of large scale (Moffat 1994). Such rare large clumps/shocks also emit and scatter more photons than their smaller cousins, which may lead to detected stochastic fluctuations, as now seen in ζ Pup in the optical with BRITE, and possibly in X-rays in the observations described here.

When trying to determine the physical origin of some stable monoperiodicity in the signal received from a star, rotation is of course the most obvious culprit. In addition, the fact that we have now clearly seen the same period in two wavelength bands has other ramifications. If rotation is accepted as the principal motor driving a periodic behavior of the star, there must be some way of connecting surface inhomogeneties in the optical photosphere with some sort of modulation of the X-ray emitting and absorbing regions of the wind. This makes it very tempting to invoke CIRs, with the optical (continuum) light curve arising from the bright spot on the stellar surface at the base of the CIR and the X-ray light curve arising from time-varying visibility of shocks in the CIR due to its rotational modulation out in the wind. Qualitatively, this idea fits with intermediate phase lags up to 0.1 in phase seen in three different successive optical recombination lines compared to the continuum light coming from the photosphere (Ramiaramanantsoa et al. 2018). There are two possible mechanisms that could cause a CIR imprint on the X-ray signal: (1) as a CIR sweeps through the unocculted portions of the wind, rotation could modulate the X-ray signal by revealing zones of additional emission, or (2) obscuration of the shocks could be caused by increasingly more or less wind material between the X-ray sources and the observer. If CIRs are indeed causing the modulation of the X-ray signals, there is an interesting interplay between the inclination angle of the star and the physical structure of CIRs. Obviously, if the star is seen exactly pole-on, there would be no rotational modulation of the signal. If the inclination angle of the star could be independently constrained, that knowledge would reveal information about the extent of CIRs in both radius and latitude. At small inclination angles, only X-ray-emitting material very near the star would be noticeably occulted each rotation.

In addition to the X-ray modulation described here, DACs, thought to be an observational manifestation of CIRs, have been observed in the past for this star in the UV. This independent line of evidence indicates that the inclination of ζ Pup must be far enough from pole-on to allow some rotational modulation of the signal. If a low inclination means our line of sight to this star may just clip the high-latitude portions of the CIRs, this may explain why DACs have been more variable in this star than in some other similar stars.

It should be noted that while it is difficult to explain any regular periodicity with periods longer than the rotational period, it is relatively easy to explain any shorter period harmonically related to the rotation period of the star. Physically, this would be manifested as having multiple (more or less persistent) structures at different longitudes around the star. If, for instance, there are two CIR footpoints ("hot spots"), the rotation period P_rot could be 1.78 day or 2 × 1.78 days = 3.56 days. This configuration (two structures per 2π) has been discussed as being plausible by Kaper et al. (1999), de Jong et al. (2001), and Massa et al. (2019), and is consistent with the period we determined of 1.7820 days. However if the period is 3.56 days, the two structures must be sufficiently symmetrical that the true rotational fundamental of 3.65 days does not show up strongly in our data or in optical data. Further knowledge of this star's CIRs is needed to gauge the likelihood of this scenario of very symmetric CIRs separated by exactly 180°. Two rotation possibilities related to the 1.78 day signal can be summarized as follows:

• 1.  
  If P_rot=1.78 days, the star would rotate very close to breakup, which poses questions regarding the wind driving, the wind symmetry, and the lack of a disk. A star so near breakup would be somewhat oblate. One or more hot spots could be present.
• 2.  
  If P_rot  = 2 × 1.78 day, the question is why this period does not show up in the periodogram above noise level. A significant degree of symmetry would be needed between the two hemispheres. This scenario would imply rotation velocity safely below critical and perhaps a more equator-on view. While this rotation period would still comply with the limits posed by $v\sin i$ and distance, the fit is not as comfortable as for the 1.78 day period.

It might be difficult to discriminate between these two cases because both of these possible rotation periods fit (if just barely for the longer one) in the range of allowable rotation periods in a recent detailed analysis of the fundamental properties of this star. Howarth & van Leeuwen (2019) argue that the "revised" Hipparcos distance of d = 332 ± 11 pc is reliable. Using this as their basis, their analysis went on to exclude any period above 3.7 days with 95% confidence. The argument was based on the robust spectroscopic measurements of $v\_\mathrm{eq}\times \sin (i)=213\,\pm 7$ km s⁻¹, and the obvious maximum of $\sin (i)=1$ (equator-on view).

Howarth & van Leeuwen (2019) discuss the consequences of the surprisingly small Hipparcos distance on ζ Pup physical properties, namely that it would have a relatively low luminosity and thus mass, atypical for its spectral type. If ζ Pup in fact has more normal properties for its spectral type, its true distance is larger and a rotation period of 2 × 1.78 day is plausible within the limit posed by the measured $v\sin i$.

If we relax some of these assumed constraints, it is within the realm of possibility that the 5 day or 6 day period apparent in our data is the true rotational period of this star. We however consider this unlikely for the following reasons. First, to have such a long rotational period, the distance to the star must be very much greater than that found by Howarth & van Leeuwen (2019) as described above. Second, it is not clear why the 1.78 day signal would be so strong when it would just be a harmonic of the fundamental (rotation) period. To hypothesize a rotation period in the 5 day or 6 day range, the 1.78 day period would need to be the n = 3 harmonic of the fundamental, indicating that the true rotation period would be 5.35 days. The possibility of the 1.78 day period being the n = 3 harmonic is discussed in Section 3.2. If the 1.78 day period were the 4th harmonic of the rotational period, the rotational period would have to be 7.12 days, outside the range of periods found on the periodogram, and this is of course strongly excluded under the recent distance determination.

To combine these effects into a specific physical example, let us say that the rotational period of the star is indeed three times the 1.78 day period (5.35 days). This would imply a distance of 480 pc or more, which does not look likely in light of the recent work of Howarth & van Leeuwen (2019) described above. We therefore conclude that the "rotation period" (as defined above) is most likely either 1 × 1.78 days or 2 × 1.78 days.

Finally, there is a note of caution that should be applied when comparing periods (and harmonics) measured for an individual star using different methods and at different epochs. There is no a priori reason why hot stars should be solid body rotators. Howarth & van Leeuwen (2019) proposed anti-solar differential rotation. Structures at different latitudes could be going around at slightly different angular speeds. Alternatively, a structure causing variability in some specific wave band could migrate in latitude over time, causing a change in the measured rotation period. In our case, P_rot means the rotation period at the latitude where the structures are located that give rise to the observed variability at the epoch of observation.

When evaluating possible physical explanations for these patterns of variability, it is useful to review the evolutionary history of ζ Pup. ζ Pup is currently believed to be a single, massive runaway star. When it was in a binary system with a more massive primary star, the Roche Lobe Overflow (RLOF) process would have spun up ζ Pup in the time before the primary's SN explosion. After that explosion, the secondary star (now ζ Pup) was ejected from the system with high spin rate in the opposite direction to that of the remnant compact primary. It therefore would not be surprising to find that ζ Pup shows rapid rotation. It was noted above that accepting a relatively close distance for ζ Pup would indicate that it is under-luminous for a star of its spectral type, but it is possible to explain this discrepancy by appealing to its individual evolutionary history in a mass-exchange binary.

The putative evolutionary history of ζ Pup could conceivably make a different contribution to the periodic variability of this system. There is some finite possibility that ζ Pup remained bound with at least a part of the debris from the exploding star or material participating in the RLOF. The 5 day or 6 day periodicity could be caused by some sort of low-mass companion object (of whatever origin) orbiting ζ Pup at a distance of ∼3 stellar radii. Rotation is the obvious prime mover for any clock-like periodicity for this star, but orbital motion provides many other options. Orbital motion could explain any stable periods longer than the rotational period. Such orbital motion periods would be expected to have no specific relation to that rotational period. While this scenario is highly speculative, we mention it for the sake of completeness.

An additional source of periodicity could be Nonradial Pulsations (NRPs). As discussed in Howarth & van Leeuwen (2019), Howarth & Stevens (2014) had applied the theoretical pulsation models of Saio (2011) to this star. Acceptance of the new nearer distance, with its attendant under-luminosity, and the individual evolutionary history of ζ Pup make it difficult to apply the Saio (2011) models, which are for standard single-star evolutionary tracks. Without detailed modeling based on individual properties of ζ Pup giving specific modes, periods, and amplitudes, it is not possible to evaluate what contributions these pulsations make to the observed periodic variabilities.

Setting aside the determination of the causes of the specific values of the periodicities in the observation, the other extremely interesting aspect of the data is the phase difference between the flux maxima in the X-ray and optical folded light curves. This paper will not attempt detailed modeling, but from a theory perspective a most important clue to the nature of the variability is the fact that most of the roughly 6% (peak-to-peak) X-ray variation is coherent on a period of 1.7806 days over a time of up to a year. If we assume that the maximum of the BRITE light curve at ϕ = 0.45 is related to the Chandra light-curve maximum at ϕ = 0.9, the time lag is about 6.1 × 10⁴ s relative to the BRITE maximum. If the 1.78 day period is interpreted as a rotation period (as in Ramiaramanantsoa et al. 2018), this corresponds to a phase lag of some 45% of a cycle between the optical maximum and the X-ray peak. A potential interpretation for this phase lag is the curved shape of a CIR in the wind (Cranmer & Owocki 1996), caused by the effects of rotation on wind streams with different rates of radial acceleration.

Analyzing the coriolis influence on wind acceleration in the corotating frame shows that the radius of a CIR shock, in stellar units, caused by a bright photospheric spot (Cranmer & Owocki 1996), is characterized by the ratio of the terminal speed to the rotation speed, times the angular scale (in radians) of the spot on the stellar surface. The phase lag of the X-ray hot spot where streamlines converge is of the same order as the size of the spot. The spot brightness contrast must be fairly significant in order to produce a strong shock (by overloading the local mass flux, it stalls and gets rammed from below). Because the optical brightness variations are at only the 1% level, it means the spot must be relatively small, no more than a few percent of the stellar surface and covering a phase interval no larger than 0.1 of the rotation period. For ζ Pup the ratio of terminal speed to rotation speed is about 10, so the shock forms at much less than 10 stellar radii, possibly even in the range of 1–3, where we also expect the bulk of the X-rays to form. Because small spots and low X-ray radii are consistent with a lag in the X-ray emission of no more than about 0.1 of the rotational period, we would have difficulty forming a consistent picture if we thought the X-ray peak was associated with the BRITE peak about ϕ = 0.45 earlier in the rotational period. Thus, either the X-ray peak is associated with the weaker BRITE peak near phase 0.8, or the BRITE signal is created by a dark spot offset by ϕ = 0.5 from our expectation. Although dark spots are equally capable of producing X-ray shocks, because the latter require only a change in terminal speed and dark spots can underload the wind and generate fast streams, Cranmer & Owocki (1996) found that dark spots do not generate DACs. Hence, the most self-consistent interpretation is that the X-ray peak is associated with the second, albeit weaker, BRITE peak. What would be necessary to verify this interpretation is simultaneous X-ray, optical, and UV line observations, to test the rotational phase relationships of DACs and CIRs relative to X-ray generation.

Thus, it seems more likely that the two optical light-curve maxima represent two different surface hot spots on the star. The Chandra light curve also has two potential maxima. The similarity between the phase lags for each of these X-ray/optical pairs suggests that, rather than a phase difference of ϕ = 0.45 for the primary maxima only, there are two optical features, each with an associated X-ray peak and a lag time of ∼ ϕ = 0.1. This value is similar to lag times for other stars comparing X-ray and multiwavelength light curves (Massa et al. 2019). If this is in fact the case, the situation is that the largest optical maximum is associated with a rather small secondary maximum in X-rays, while the secondary optical maximum is associated with the largest feature in the X-ray light curve. We can perhaps understand such a configuration if we consider that the X-ray flux, as mentioned above, depends on the viewing angle of the curved CIR in the wind and on occultation. Thus, the strength of the X-ray signal may not be clearly correlated to the structure of the optical emission from the hot spot itself. A DAC period of about 20 hr could be explained if the light curves are interpreted as displaying evidence for two hot spots per rotation cycle of ∼2 × 20h = 40h. DAC periods might fit into this scenario with periods of about 0.8 days, which is about half of the rotation period.