The distance, rotation, and physical parameters of zeta Pup

TL;DR

This paper confirms the distance to zeta Pup, refines its physical parameters, and discusses implications for its rotation, evolution, and wind variability, suggesting a binary interaction history and challenging some models of stellar wind behavior.

Contribution

It provides a reliable distance measurement, constrains the star's rotation period, and links physical properties to binary evolution and wind variability models.

Findings

Distance to zeta Pup is 332+/-11 pc.

Rotation period is less than 3.7 days, ruling out ~5.1 days.

Physical parameters are below canonical values, indicating binary interaction history.

Abstract

We scrutinize the Hipparcos parallax for the bright O supergiant zeta Pup, and confirm that the implied distance of 332+/-11pc appears to be reliable. We then review the implications for the star's physical parameters, and the consequences for the interpretation of P(phot), the 1.78-d photometric period. The equatorial rotation period is <3.7d (with 95% confidence), ruling out a proposed ~5.1d value. If the photometric period is the rotation period then i, the inclination of the rotation axis to the line of sight, is 33.2+/-1.8 degrees. The inferred mass, radius, and luminosity are securely established to be less than canonical values for the spectral type, and are not in agreement with single-star evolution models. The runaway status, rapid rotation, and anomalous physical properties are all indicative of an evolutionary history involving binary (or multiple-star) interaction. We perform simple starspot modelling to show that the low axial inclination required if P(rot) = 1.78d has testable spectroscopic consequences, which have not been identified in existing time series. If P(phot) is directly related to drivers of systematic, high-velocity stellar-wind variability (`discrete absorption components') in zeta Pup, antisolar differential rotation is required. Model line profiles calculated on that basis are at variance with observations.

Figures

Tables

Table 1. Table 1: Selected model results. The second column summarizes the discussion of Section  4.1 . Subsequent columns report the more through analyses described in Section  4.2 , under the assumptions stated in the column headers. The equatorial Eddington gravity, g e Edd subscript superscript 𝑔 Edd e g^{\rm Edd}_{\rm e} , is defined in Section  4.4 .

	Spherical			$i={90}^{\circ }$			$P_{\mathrm{rot}}=1.78$ d
$M/{\text{M}}_{\odot }$	25.3	$\pm$	5.3	26.6	$\pm$	5.6	$\equiv 15$			$\equiv 25$			$\equiv 50$
$R_p$/${\text{R}}_{\odot }$	13.50	$\pm$	0.52	13.22	$\pm$	0.54	10.06	$\pm$	0.07	11.25	$\pm$	0.19	12.32	$\pm$	0.34
$R_e$/${\text{R}}_{\odot }$				14.08	$\pm$	0.53	13.86	$\pm$	0.05	13.72	$\pm$	0.49	13.66	$\pm$	0.54
$T_p/$ kK‡	40			41.59	$\pm$	0.32	46.21	$\pm$	0.38	44.35	$\pm$	0.34	42.47	$\pm$	0.25
$T_p/T_e$	$\equiv 1.0$			1.064	$\pm$	0.014	1.392	$\pm$	0.064	1.207	$\pm$	0.023	1.103	$\pm$	0.012
${\log }_{10}{(L/{\text{L}}_{\odot })}^{‡}$	5.65	$\pm$	0.06	5.641	$\pm$	0.033	5.536	$\pm$	0.019	5.575	$\pm$	0.024	5.602	$\pm$	0.030
${\log }_{10}\left(g_p\right)$, cgs	3.58	$\pm$	0.08	3.611	$\pm$	0.078	3.609	$\pm$	0.006	3.734	$\pm$	0.014	3.956	$\pm$	0.024
${\log }_{10}\left(g_e\right)$, cgs				3.494	$\pm$	0.103	2.695	$\pm$	0.189	3.309	$\pm$	0.069	3.761	$\pm$	0.048
${\log }_{10}\left(g_e^{\mathrm{Edd}}\right)$, cgs	3.23			3.19			2.905	$\pm$	0.065	3.080	$\pm$	0.019	3.159	$\pm$	0.008
${\omega }_e/{\omega }_{\mathrm{crit}}$	—			0.602	$\pm$	0.053	0.985	$\pm$	0.010	0.902	$\pm$	0.022	0.731	$\pm$	0.030
Inclination $i{(}^{\circ })$	—			$\equiv 90$			32.8	$\pm$	1.7	33.2	$\pm$	1.8	33.4	$\pm$	1.9
$v_e$ (km s-1)	213	$\pm$	7	213	$\pm$	7	394	$\pm$	14	390	$\pm$	14	388	$\pm$	15
$P_{\mathrm{rot}}$ (d)	$\le 3.21$	$\pm$	0.17	3.35	$\pm$	0.16	$\equiv 1.78$
‡Assuming $T_{\mathrm{eff}}=40$ kK, excepting the spherical-star luminosity; if $f=T_{\mathrm{eff}}^{\ast }/40$ kK, where $T_{\mathrm{eff}}^{\ast }$ is the true effective temperature, then $T_p^{\ast }=f{T}_p$ and ${\log }_{10}\left(L^{\ast }/{\text{L}}_{\odot }\right)={\log }_{10}\left(L/{\text{L}}_{\odot }\right)+4{\log }_{10}f$. Radii at $P_{\mathrm{rot}}=1.78$ d are smaller than for the $i={90}^{\circ }$ models for the reason mentioned in Section 4.3 – the hotter polar regions are more clearly presented to the observer at $i\simeq 33$ than equator-on, so a smaller emitting area is required to match the observed flux, for given $T_{\mathrm{eff}}$.

Table 2. Table 2: Bonnsai model results (cf. Section  4.5 ). Small open circles indicate parameters used as inputs (with observational values and 1- σ 𝜎 \sigma gaussian errors listed in column 2); quantities assessed on the ZAMS are indicated with a subscript ‘1’.

Parameter	Input			Replicated Observables
				Model 1			Model 2
$T_{\mathrm{eff}}$/kK	40	$\pm$	1	$\circ$	$39.92$	${}_{-0.88}{}^{}{}_{}{}^{+1.16}$	$\circ$	$40.08$
${\log }_{10}\left(L/{\text{L}}_{\odot }\right)$	5.60	$\pm$	0.05	$\circ$	$5.58$	${}_{-0.05}{}^{}{}_{}{}^{+0.05}$	$\circ$	$5.58$	${}_{-0.04}{}^{}{}_{}{}^{+0.06}$
$v_e\sin i$/km s-1	213	$\pm$	7	$\circ$	$210$	${}_{-6}{}^{}{}_{}{}^{+12}$	$\circ$	$210$	${}_{-5}{}^{}{}_{}{}^{+13}$
$Y$	0.41	$\pm$	0.05		$0.27$	${}_{-0.01}{}^{}{}_{}{}^{+0.04}$	$\circ$	$0.39$	${}_{-0.06}{}^{}{}_{}{}^{+0.09}$
$v_e$/km s-1					$220$	${}_{-15}{}^{}{}_{}{}^{+58}$		$270$	${}_{-17}{}^{}{}_{}{}^{+36}$
$v_{e,1}$/km s-1					$260$	${}_{-29}{}^{}{}_{}{}^{+114}$		$440$	${}_{-40}{}^{}{}_{}{}^{+47}$
$M_{\mathrm{evol}}/{\text{M}}_{\odot }$					$42.4$	${}_{-3.6}{}^{}{}_{}{}^{+3.4}$		$37.0$	${}_{-2.9}{}^{}{}_{}{}^{+3.8}$
$M_1/{\text{M}}_{\odot }$					$45.4$	${}_{-3.6}{}^{}{}_{}{}^{+3.1}$		$41.8$	${}_{-2.7}{}^{}{}_{}{}^{+3.6}$
Age/Myr					$2.20$	${}_{-0.49}{}^{}{}_{}{}^{+0.57}$		$3.56$	${}_{-0.75}{}^{}{}_{}{}^{+0.77}$

Equations

$$f_V=(5.14\pm 0.17)$$

$$[\equiv (5.14\pm 0.17)$$

$$F_V(=4\pi H_V)=(6.13\pm 0.13)\times 10^8\text{ erg cm−2 s−1 A˚−1}.$$

$$R_{\mathrm{eff}}=13.50\pm 0.52{\text{R}}_{\odot };$$

$$\log (L_{\mathrm{eff}}/{\text{L}}_{\odot })=5.65\pm 0.06$$

$$P_{\mathrm{rot}}^{\max }=2\pi R_{\mathrm{eff}}/v_{\mathrm{e}}\sin i=3.21\pm 0.17\text{ d},$$

$$\sin (i)=(P_{\mathrm{rot}}{v}_{\mathrm{e}}\sin i)/(2\pi R_{\mathrm{eff}}).$$

$$M=22.1\pm 4.6{\text{M}}_{\odot }.$$

$$g_{\mathrm{N}}\simeq g_{\mathrm{obs}}+(v_{\mathrm{e}}\sin i{)}^2/R_{\ast }$$

$$M=25.3\pm 5.3{\text{M}}_{\odot }.$$

$$Q$$

$$=0.180\pm 0.018\text{ d}$$

$${\omega }_{\mathrm{crit}}=\sqrt{(GM)/(1.5R_{\mathrm{p}}{)}^3}$$

$$M\simeq g_{\mathrm{N}}{R}_{\mathrm{p}}{R}_{\mathrm{e}}/G$$

$$T_{\mathrm{eff}}=\sqrt[4]{\int \sigma (T_{\mathrm{eff}}^{\ell }{)}^4\text{d}A/\int \sigma \text{d}A}$$

$$g_{\mathrm{Edd}}$$

$$\simeq 6.56\times 10^{-16}(T/\text{K}{)}^4{\text{ cm s}}^{-2}$$

$$\frac{\omega (\theta )}{{\omega }_{\mathrm{e}}}$$

Full text

††pubyear: 2018††pagerange: The distance, rotation, and physical parameters of $\zeta$ Pup–References

1 Introduction

As the brightest early-O supergiant by a comfortable margin, $\zeta$ Pup (HD 66811, O4 I (n)fp; Sota et al. 2014) has long been a touchstone in the development of models of radiatively-driven winds, and their role in massive-star evolution. In parallel, it has been subject to observational scrutiny at all accessible wavelengths, from X-ray to radio (e.g., Blomme et al. 2003; Eversberg et al. 1998; Hanson et al. 2005; Harries & Howarth 1996; Marcolino et al. 2017; Nazé et al. 2018; Reid & Howarth 1996; and many others). However, it was only with the relatively recent availability of satellite photometry that an apparently periodic signal with $P_{\rm phot}=1.78$ d was discovered in its optical brightness (Howarth & Stevens, 2014; Ramiaramanantsoa et al., 2018).

The origin of this signal remains moot. In a major study built on photometry obtained with BRITE-Constellation nanosatellites (Weiss et al., 2014), Ramiaramanantsoa et al. (2018) argued for it to be the rotation period, an interpretation that Howarth & Stevens (2014) had considered less likely than a pulsational origin, on the grounds that exceptional, near-critical rotation would be implied – the projected equatorial rotation velocity of $\zeta$ Pup, $v_{\rm e}\sin{i}\simeq 213$ km s*-1* ($\S$4.1.4), is already the most rapid known for any Galactic O supergiant (Howarth et al., 1997), and a 1.78-d rotation period would require that $v_{\rm e}$ be $\sim{2}\times$ greater. Nevertheless, this is merely a plausibility argument, and exceptionally rapid rotation could simply be a signature of an exceptional evolutionary history (cf., e.g., Vanbeveren 2012; de Mink et al. 2013).

The case of $\zeta$ Pup therefore represents a modern manifestation of the long-standing difficulty in making a compelling distinction, observationally, between pulsational and rotation modulation as the mechanism responsible for low-amplitude spectroscopic and photometric variability in early-type stars (cp., e.g., Gies 1991; Harmanec 1999). This problem has been compounded in the case of $\zeta$ Pup by some contention in respect of its distance, and consequently in basic physical parameters such as radius and mass. Our purpose in this paper is to examine these issues.

To do so, we first review the Hipparcos parallax data in Section 2 (and examine the Hipparcos photometry in Section 3). We evaluate basic stellar parameters in Section 4, under the limiting assumptions of (i) axial inclination $i=90^{\circ}$ and (ii) rotation period $P_{\rm rot}=P_{\rm phot}(=1.78$ d). Evolutionary implications are discussed in Section 4.5. Physical modelling of photometric and spectroscopic variability, intended to test the rotational-modulation/hotspots hypothesis, is presented in Section 5. The case for an association between $P_{\rm phot}$ and discrete absorption components in the stellar wind is scrutinized in Section 6.1. Finally, the summary and conclusions are given in Section 7.

2 Distance

2.1 Hipparcos

Although the revised reduction of the Hipparcos data yields a reasonably precise parallax of $(3.01\pm 0.10)$ mas for $\zeta$ Pup (distance $d=332\pm 11$ pc; van Leeuwen 2007a, b), subsequent state-of-the-art model-atmosphere analyses have disregarded or challenged this result (which implies unexpectedly low values for the stellar mass and radius,111Equivalently, the implied absolute magnitude, $M(V)=-5.5$ ($\S$4.1), is $\sim$1m fainter than canonical values for the spectral type (e.g., Walborn 1973; Martins & Plez 2006), and is close to expectations for a main-sequence star. for standard evolutionary scenarios), preferring larger distances (up to $\sim$700 pc; e.g., Najarro et al. 2011; Bouret et al. 2012; Pauldrach et al. 2012), with concomitant implications for the luminosity, etc.

Unfortunately, at $V\simeq 2.1$ $\zeta$ Pup is too bright to have been included in currently available Gaia data releases. However, its brightness is an asset where Hipparcos is concerned (parallax errors are limited by photon noise at $V\sim 3$ and fainter, although by calibration uncertainties otherwise222Compared to the original analysis, calibration uncertainties were reduced by a factor $\sim$5 in the 2007 re-reduction.), and we have reviewed the results to check if there are any reasons to suspect the published parallax.

We find no suggestion of any problems in the astrometric data; the error correlations for the astrometric parameters are very low, and the distribution over scan directions very good, as is the distribution over parallax factor. The underlying data for the star are consistent and numerous (138 observations, with only 4 rejections in the iterative solution); and $\zeta$ Pup is in a part of the sky where the scan coverage is almost maximally good.

[The parallax factor multiplied by the actual parallax of the star gives the along-scan displacement of the position due to the parallax at the time of observation, for the given scan direction; it always lies in the range $-0.7$:+0.7. For $\zeta$ Pup the parallax factors fall entirely in the ranges $-0.7$:$-0.4$ and +0.4:+0.7, which is very good distribution for a reliable parallax determination. The spread over epochs is also very good, giving a low ($\sim\mathcal{O}$(0.1)) error correlation between proper-motion and parallax determinations.]

As a further check, we compared results for apparently single stars within 180′ of $\zeta$ Pup that have measurements both from Hipparcos and in Gaia DR2 (Gaia Collaboration, 2018a). Fig. 1 shows the result. This confirms that there are no reasons to doubt the Hipparcos measurement.

Finally, given that the parallactic displacement is $\sim$10$\times$ greater than the star’s angular diameter (Section 4.1.3), and that the amplitude of photometric variability is small, any asymmetry in the surface-brightness distribution is very unlikely to compromise this conclusion.

2.2 Corroboration

Since the discovery of the Gum nebula, its ionization has generally been attributed to $\gamma^{2}$ Vel and $\zeta$ Pup (Gum, 1952). Prior to the advent of satellite astrometry, the distance to $\zeta$ Pup was therefore estimated on that basis (e.g., “both $\gamma^{2}$ Vel and $\zeta$ Pup appear to be embedded in the giant Gum H ii region and are the sources of its ionization so that we can assume the two stars are at the same distance”; Morton et al. 1969).

In support of this assumption, Woermann et al. (2001) argued that the Gum nebula could be the remnant of a supernova from a previous binary companion to $\zeta$ Pup, noting that the surviving O star passed within $\nicefrac{{1}}{{2}}^{\circ}$ of the expansion centre of the nebula about 1.5 Myr ago. They further concluded that the nebula is probably primarily ionized by $\zeta$ Pup, at a distance in the range $\sim$200–500 pc.

A physical association of the various components of the ‘Vela Complex’, including the Vela supernova remnant, $\gamma^{2}$ Vel and $\zeta$ Pup, the Gum nebula, and the Vela OB2 association, has also been widely assumed (e.g., Sushch et al. 2011). The significance of this is that apparently reliable distant estimates for other components of the Vela Complex can be used as a check on the plausibility of the Hipparcos parallax for $\zeta$ Pup.

The $\gamma^{2}$ Vel binary system is particularly useful in this context, as its distance can be independently established by primary (geometric) means. Millour et al. (2007) obtained a distance $d=368^{+38}_{-13}$ pc by combining new interferometric separation measurements with spectroscopic-orbit observations, a result independently confirmed and refined by North et al. (2007; cf. also Lamberts et al. 2017), who, in effect, solved the orbit in three dimensions to obtain $d=336^{+8}_{-7}$ pc. The new reduction of Hipparcos data for $\gamma^{2}$ Vel gives a distance $d=342^{+40}_{-32}$ pc (van Leeuwen, 2007a, b), in excellent agreement. Further corroboration is provided by a precise photometric determination of the distance to Vela OB2, yielding $d=350\pm 13$ pc (Jeffries et al., 2009).

We conclude that independent determinations of the distances to Vela OB2 and to $\gamma^{2}$ Vel are in good mutual agreement, and both are in good accord with the Hipparcos distance to $\zeta$ Pup. Once again, there appear to be no good grounds to doubt the reliability of the Hipparcos parallax for $\zeta$ Pup.

2.3 The runaway $\zeta$ Pup

That $\zeta$ Pup is a runaway star was first proposed by Upton (1971), and its dynamics and origin have subsequently been discussed a number of times (e.g., Blaauw 1993; van Rensbergen et al. 1996; Moffat et al. 1998; Hoogerwerf et al. 2001; Schilbach & Röser 2008).

The Hipparcos  parallax and proper motion ($34.07\pm 0.10$ mas yr*-1*) yield a transverse velocity of $53.7\pm 1.8$ km s*-1*. We found four primary literature sources that yield 19 separate radial-velocity measurements (Frost et al. 1926; Wilson 1963; Conti et al. 1977; Garmany et al. 1980); those measurements are in satisfactory mutual agreement, with roughly similar estimated errors (and provide no evidence for binarity). We adopt their unweighted average, $-25.4\pm 2.1$ km s*-1* (s.e.).

The space velocity with respect to the Sun is thus $59.4\pm 1.9$ kms. Correcting for the Sun’s peculiar motion and for Galactic rotation333Using (U,V,W){}_{\odot}=(11.10,12.24,7.25)~{}\mbox{\rm km;$s${}^{-1}} (Schönrich et al., 2010); $R_{0}=8.5$ kpc, \theta_{0}=220~{}\mbox{\rm km;$s${}^{-1}} (Kerr & Lynden-Bell, 1986); and the Galactic pole position adopted in ESA (1997). Quoted errors do not include uncertainties on these quantities. we obtain a peculiar velocity of 56.2\pm 1.9~{}\mbox{\rm km;$s${}^{-1}} (transverse and radial components 36.6, $-42.7$ km s*-1*). These figures confirm that $\zeta$ Pup is a runaway by any generally accepted definition (e.g., Gies 1987).

However, the Hipparcos distance rules out the runaway scenarios discussed by van Rensbergen et al. (as already noted by Schilbach & Röser 2008), and hence also the specific binary-merger evolutionary scenario proposed by Vanbeveren (2012) and discussed by Pauldrach et al. (2012). The most likely site of origin for $\zeta$ Pup appears to be the cluster Trumpler 10 (Hoogerwerf et al., 2001; Schilbach & Röser, 2008) on the basis of its Hipparcos distance ($d=386\pm 5$ pc; van Leeuwen 2009), although the Gaia value ($441\pm 4$ pc; Gaia Collaboration 2018b), based on a much larger sample of stars, may require a review of this conclusion.

3 Hipparcos photometry

As well as astrometry, Hipparcos provided relatively precise broad-band photometry; results for $\zeta$ Pup are shown in Fig. 2. The formal errors may not be reliable for such a bright target, due to saturation effects, but nevertheless the dispersion of the data (s.d. 7.5 mmag) is consistent with the stochastic microvariability on $\sim$10-hr timescales reported by Ramiaramanantsoa et al. (2018; see also Balona 1992), with possible longer-term changes at the $\sim$1% level.

These data were previously examined by Marchenko et al. (1998), who found a 2.563-d periodicity with semiamplitude 6 mmag (again, tentatively attributed to rotational modulation). Our time-series analysis of these data recovers this result, but shows no evidence for a significant periodic signal close to $P_{\rm phot}=1.78$ d. Fig. 2 includes the photometry phase-folded on an ephemeris with $T_{0}=\text{JD }2\ 448\ 464.0$ (close to the mean date of observation), $P=1.780\ 938$ d (Howarth & Stevens, 2014), and binned in phase intervals of 0.1. Error bars shown on the binned points are standard errors, using standard deviations in each bin computed from the dispersions of the data (and not the nominal errors on individual points).

These binned data are not formally consistent with a phase-constant flux value, with $\chi^{2}=38.7$, but this rather large value is dominated by a single outlying bin, at $\phi\simeq 0.45$ in Fig. 2. Of only five observations falling in that bin, four accordant results come from a single epoch; the 5th observation is well separated in time, but differs by less than 1 mmag from the mean brightness, resulting in a very small error bar. We conclude that the phased data are consistent with a near-constant time-averaged flux (semi-amplitude $\lesssim$3 mmag) over at least 80% of the 1.78-d period.

A 1.78-d periodic signal in the photometry could be smeared out by phase drift. However, examination of subsets of the epoch photometry shows no evidence for short-lived, large-amplitude periodic signals. Moreover, the SMEI discovery data have a very similar time-span to the Hipparcos mission ($\sim$1000 d), yet show obvious, near-coherent, $\sim$sinusoidal changes (Fig. 2). Additionally, Balona (1992) failed to detect a signal in precise ground-based photometry obtained in 1989 April, $\sim$four months before the first Hipparcos observations. A plausible interpretation is that, around $\sim$1990, the 1.78-d signal was of significantly lower amplitude than in $\sim$2005 and 2015, when a $\gtrsim$5 mmag semi-amplitude was recorded (Howarth & Stevens, 2014; Ramiaramanantsoa et al., 2018).

4 Rotation and physical parameters

Given the apparently well-determined distance, we can review the implied physical parameters for $\zeta$ Pup, including its rotation (with the ultimate aim of testing the hypothesis that the 1.78-d photometric signal is a direct tracer of the rotation period).

4.1 Principles & preliminaries: slow-rotation limit

The basic principles are simple: the observed and surface fluxes lead to an estimate of the angular diameter; the angular diameter and distance give the radius; the radius and an assumed rotation period lead to the equatorial rotation speed, $v_{\rm e}$; which, with the observed projected rotation speed, gives the inclination of the rotation axis to the line of sight, $i$. The radius, effective temperature, and surface gravity also yield the luminosity and (‘spectroscopic’) mass.

These are straightforward sums as long as the star is satisfactorily approximated as a sphere of uniform surface flux. This is marginal for $\zeta$ Pup; its rapid rotation introduces complications that are considered further in Section 4.2. Nevertheless, an examination of parameters in the spherical-star, or slow-rotation, limit is of use to provide a context and a point of comparison with previous analyses, and allows us to assemble some necessary numerical data.

4.1.1 Observed flux

The observed $V$ magnitude and colours yield the extinction. We took $V=2.25$, $(B-J)=-0.82$ (Johnson et al., 1966) to estimate $E(B-V)=0.04$, and thence a reddening-corrected visual magnitude $V_{0}=2.13$. Using the Hipparcos parallax, the implied absolute magnitude is $M_{V}=-5.48\pm 0.09$.

The corresponding flux is

[TABLE]

Here we have used an extinction law, intrinsic colours, and flux calibration from Howarth (1983, 2011). The error quoted on the observed flux is based on 1% uncertainties in each of $V$, $E(B-V)$, and the absolute calibration.

4.1.2 Surface flux

A number of nLTE model-atmosphere analyses of $\zeta$ Pup have been published; we use results from Kudritzki et al. (1983), Bohannan et al. (1986), Puls et al. (2006, supplanting Repolust et al. 2004), Bouret et al. (2012), and Pauldrach et al. (2012), each of whom employed different, independent modelling codes.444The two 20th-century analyses are based on plane-parallel, hydrostatic models, without line blanketing; their surface fluxes bracket those of the 21st-century studies, which allow for stellar winds and line blanketing. Each study also made different assumptions about the distance to $\zeta$ Pup, but we can take the authors’ adopted distances and derived radii,555In every case, the authors assumed spherical symmetry. together with the observed flux, to infer their model’s $V$-band surface flux. We find

[TABLE]

Here, as in most of the subsequent analysis, we adopt the standard deviation of model-atmosphere results as a more conservative, and arguably more credible, estimate of the true uncertainty associated with the spectroscopic analyses than is provided by the formal standard error. In this case the s.d.  reflects differences in different modelling procedures, input physics, and numerical methods, as well as differences in inferred atmospheric parameters (principally $T_{\rm eff}$, but also $\log_{10}(g)$, helium abundance, etc.).

4.1.3 Effective radius, luminosity

The observed and model fluxes yield the effective angular diameter directly: $\theta_{\rm eff}=2\sqrt{f/F}=0.378\pm 0.007$ mas, consistent with the observed value of $0.41\pm 0.03$ mas666Reduced from the published limb-darkened value of 0.42 to correct for electron scattering in the wind (Kudritzki et al., 1983). We note, however, that formally statistically significant discrepancies between the pioneering intensity-interferometer results and modern long-baseline optical interferometry are not uncommon (cf. Baines et al. 2018). (Hanbury Brown et al. 1974; here we use the ‘effective’ qualifier to indicate the result for a spherical star, of uniform surface flux, that matches the observed $V$-band brightness of $\zeta$ Pup).

The angular diameter and parallax may be combined to give the effective radius,

[TABLE]

the error estimate takes into account uncertainties in the parallax, the observed reddening-free flux, and the surface flux. Rescaling the luminosities from the nLTE analyses listed in Section 4.1.2 to the Hipparcos distance we find a corresponding effective luminosity

[TABLE]

(where the quoted error reflects the dispersion in the analyses and the uncertainty in the distance).

4.1.4 Preliminary rotation, inclination

For a spherical star the rotation period is $P_{\rm rot}=2\pi R_{\rm eff}/v_{\rm e}$, so an upper limit, $P_{\rm rot}^{\rm max}$, follows from the radius and observed $v_{\rm e}\sin{i}$ by assuming $\sin{i}=1$.

Observational determinations of $v_{\rm e}\sin{i}$ are in remarkably good accord; the seven values independently determined by Kudritzki et al. (1983), Bohannan et al. (1986), Penny (1996), Howarth et al. (1997), Repolust et al. (2004), Bouret et al. (2012), and Pauldrach et al. (2012) all lie in the range 203–220 km s*-1*, averaging $213\pm 7$ km s*-1* (s.d).

[In principle, gravity darkening could result in all the empirical determinations of $v_{\rm e}\sin{i}$ systematically underestimating the true value (Townsend et al., 2004), but synthetic spectra from models such as those discussed in $\S$4.2 indicate that this effect is negligible for $\zeta$ Pup.]

We therefore have

[TABLE]

which securely rules out the 5.1-d rotation period proposed by Moffat & Michaud (1981).

If an assumption is made about the rotation period, then a naïve estimate of the axial inclination may be obtained instead, from the radius and observed $v_{\rm e}\sin{i}$:

[TABLE]

If $P_{\rm rot}=1.78$ d, then $\sin(i)=0.555\pm 0.028$ ($i=33{\aas@@fstack{\circ}}7\pm 1{\aas@@fstack{\circ}}9$), and $v_{\rm e}=384\pm 15$ km s*-1*, for the simple, spherical-star case.

[Alternatively, if $P_{\rm rot}=2.56$ d (the period found in Hipparcos photometry by Marchenko et al. 1998), then $\sin(i)=0.799\pm 0.040$ ($i=53{\aas@@fstack{\circ}}1\pm 3{\aas@@fstack{\circ}}9$), and $v_{\rm e}=266\pm 10$ km s*-1*.]

4.1.5 Spectroscopic mass, pulsation constant

Spectroscopic determinations of $\log_{10}(g)$ from the sources listed in Section 4.1.2 give observed values in the range 3.4–3.6 (dex cgs), averaging $3.52\pm 0.08$ (s.d.). Combining this with the effective radius established above gives a mass

[TABLE]

However, even in the spherical-star approximation, it is possible to make a statistical correction to the observed surface gravity for the effects of centrifugal forces,777Bouret et al. (2012) applied this correction to obtain their quoted $\log_{10}(g)=3.64$; we ‘uncorrected’ this, using their adopted radius, to infer an observed $\log_{10}(g)=3.61$. in order to estimate the Newtonian gravity ($GM/R^{2}$, sometimes ambiguously referred to as the ‘true’ gravity):

[TABLE]

(Herrero et al., 1992; Vacca et al., 1996; Repolust et al., 2004). For our adopted radius and $v_{\rm e}\sin{i}$ (Sections 4.1.3, 4.1.4) this leads to $\log_{10}(g_{\rm N})=3.58$, and a revised mass of

[TABLE]

If the photometric period is a pulsation period, then the pulsation ‘constant’ is

[TABLE]

[or $Q=0.260\pm 0.026$ if $P_{\rm puls}=2.56$ d].

4.1.6 Mass-loss rate

All empirical determinations of the stellar-wind mass-loss rate require information on distances; additional uncertainties arise from ignorance, in detail, of the radial density distribution (e.g., the acceleration parameter $\beta$ of a canonical velocity law, $v(r)=v_{\infty}(1-R_{*}/r)^{\beta}$; the degree of clumping, typically parameterized by a factor $f(r)=\langle\rho^{2}(r)\rangle/\langle\rho(r)\rangle^{2}$; and the overall geometry).

Puls et al. (2006) performed an extensive analysis of multiwavelength observations (H$\alpha$, IR, mm, radio), including detailed consideration of clumping; scaling their results to the Hipparcos distance yields $\dot{M}=2.6\times 10^{-6}$ $\text{M}_{\odot}$ yr*-1* ($\dot{M}\propto d^{3/2}$). The result of a methodologically independent analysis of X-ray line profiles by Cohen et al. (2010) rescales to $(2.5\pm 0.2)\times 10^{-6}$  $\text{M}_{\odot}$ yr*-1* ($\dot{M}\propto d$), in excellent accord.

4.2 Practicalities: rapid rotation

For a rotation period $\lesssim$3 d we expect rotational effects to be non-negligible. In the Roche approximation (which we adopt here, along with a default assumption of latitude-independent angular rotation velocity) both the shape of a star and the ratio of equatorial to polar gravities are determined solely by $\omega_{\rm e}/\omega_{\rm crit}$, the ratio of the equatorial angular velocity to the critical value at which the Newtonian gravitational force is matched by the centrifugal force, where

[TABLE]

for a star of mass $M$ and polar radius $R_{\rm p}$ (e.g., Collins 1963).

In the limit that $i=90^{\circ}$, then $\omega_{\rm e}/\omega_{\rm crit}\simeq 0.60$, and the rotational distortions are modest; the equator is $\sim$2.5 kK cooler than the poles, $\sim$0.12 lower in $\log_{10}(g)$, and has a $\sim$7% greater radius (Section 4.4). However, if the rotation period really is as short as 1.78 d, then $\omega_{\rm e}/\omega_{\rm crit}\gtrsim 0.90$, and the spherical-star approximation is a rather poor one, prompting a more thorough treatment.

4.2.1 Model overview

We model the rotationally distorted, gravity-darkened star by using exoBush (Howarth & Smith, 2001; Howarth, 2016). The surface geometry is that of a Roche equipotential, divided into a large number of ‘tiles’. The specific intensity (or radiance) for each tile is interpolated from a pre-computed grid of model-atmosphere results, as a function of wavelength $\lambda$, viewing angle $\mu$,888Where $\mu=\cos\theta_{\rm n}$ and $\theta_{\rm n}$ is the angle between the surface normal and the line of sight. local effective temperature $T^{\ell}_{\rm eff}$, and local effective gravity $\log_{10}(g^{\ell})$, Doppler shifted according to the line-of-sight velocity. Results for all tiles are summed, weighted by projected area, in order to generate a synthetic spectrum. The use of specific intensities means that limb darkening is taken into account implicitly, in a fully wavelength-dependent manner. Gravity darkening is modeled in the ELR formalism (Espinosa Lara & Rieutord, 2011), which gives results close to traditional von Zeipel gravity darkening (von Zeipel, 1924), but which leads to better agreement with, in particular, interferometric observations (e.g., Domiciano de Souza et al. 2014). Intensities are interpolated in the grid of hydrostatic, line-blanketed, nLTE tlusty model atmospheres described by Reeve & Howarth (2018).

4.2.2 Procedure

The calculations essentially follow the steps outlined in Section 4.1, but make full allowance for the rotational distortion and gravity darkening, using a Monte-Carlo (MC) approach to observational uncertainties.

For each MC realisation we first generate a set of values for the observed flux, surface flux, and parallax (each drawn from the corresponding mean values and errors given previously, assuming Gaussian distributions) in order to establish a value for the effective radius.

We combine this radius with a model-atmosphere \sim$$V-band surface flux999We actually use the monochromatic flux at 546.5 nm; the exact choice is of no consequence as long as a line-free wavelength is chosen, as is the case here. to generate a pseudo-‘observed’ target flux (cf. Section 4.3). A $v_{\rm e}\sin{i}$ value is also drawn from the adopted distribution.

Further steps proceed according to two alternative assumptions about the star’s rotation:

• (i)

*Assuming $i=90^{\circ}$.

*If $\zeta$ Pup is viewed close to equator-on (as has frequently been assumed, based on its exceptionally large $v_{\rm e}\sin{i}$ value), then rotational distortions are sufficiently small that it may be considered reasonable to infer a mass from the observed gravity,

[TABLE]

(where a specific ‘observed’ $\log_{10}(g_{\rm N})$ value is generated for each MC cycle, and $R_{\rm e}$ is the equatorial radius). We compute a full, rotationally-distorted, gravity-darkening model for this mass (and $i=90^{\circ}$), taking $R_{\rm p}\simeq R_{\rm eff}$ (in practice, slightly less) as a first estimate. The flux from this model will not match the target flux, because of the rotational effects now included; so we rescale the radii appropriately, recalculate the mass, and compute a new model. The process is then iterated until the model flux matches the target flux, resulting in a self-consistent pair of $R_{\rm p},M$ values that reproduce the target flux and $v_{\rm e}\sin{i}$, for the assumed inclination.

• (ii)

Assuming $P_{\rm rot}=1.78$ d.

In this case, several sequences of models are run, each characterized by a specified, assumed mass (because ‘the’ surface gravity is a poorly defined quantity for the implied rapid rotation). For each sequence, the rotation period gives an initial estimate of the inclination. As before, we generate full models, but now iterate to identify a self-consistent pair of $R_{\rm p},i$ values that reproduce the target flux and $v_{\rm e}\sin{i}$, for the assumed mass.

4.3 An ignorable aside on model-atmosphere intensities and fluxes

In principle, the ‘target flux’ used in Section (4.2.2) could be matched to the adopted value for the observed flux, simply by adjusting the model’s (global) effective temperature, defined as

[TABLE]

(where $\sigma$ is the Stefan–Boltzmann constant and the integrations are over total surface area). However, this adjustment is unnecessary as long as only the factors involving geometry are of interest (i.e., mass, radius, inclination, $\omega_{\rm e}/\omega_{\rm crit}$); the requirement then is only that a consistent temperature be adopted, not that it be ‘correct’ (beyond first order). We simply adopt $T_{\rm eff}=40$ kK, representative of results from the detailed analyses mentioned in Section 4.1.2 (even though our tlusty model fluxes will not precisely match the fluxes of the models used in those analyses).

Nevertheless, there remains a minor inconsistency in our modelling of the 1.78-d rotation constraint, which arises because the relationship between the perceived temperature and the global effective temperature varies with inclination – a gravity-darkened star will generally appear hotter if viewed pole-on than equator-on. Consequently, as the inclination changes from one iteration to the next, the perceived temperature changes, at constant $T_{\rm eff}$ (as does the computed observed flux, even at constant radius).

This could be corrected for, given an appropriate prescription for transforming between effective and perceived temperature, but one would first need to define the latter quantity (e.g., by synthesizing the full gravity-darkened spectrum and then modelling it as though it arose from a spherical star of uniform surface intensity; the result would still depend on the analysis criteria). In practice, for the rather small range of inclinations that our models are found to span, the $V$-band flux variation from this effect is negligible (acting as additional source of very-low-amplitude noise in the radius determinations).

4.4 Results

Some results of the models are summarized in Table 1. The $i=90^{\circ}$ models represent one geometrical extreme, and provide secure upper-limit values for $R_{\rm p}$, $R_{\rm e}$, and $P_{\rm rot}$, along with lower limits to $v_{\rm e}$ and $\omega_{\rm e}/\omega_{\rm crit}$. The 2-$\sigma$ upper limit on the rotation period is $P_{\rm rot}<3.7$ d, ruling out previous suggestions that the value may be $\sim$5.1 d (Moffat & Michaud, 1981).

Additional results for $P_{\rm rot}=1.78$-d models are shown in Fig. 3. These models provide some limits on the allowed masses for this rotation period. First, $M\gtrsim 12.8\text{M}_{\odot}$ is required for $\omega_{\rm e}/\omega_{\rm crit}\leq 1$. Secondly, the ‘Eddington gravity’ required to retain material against the radiation force is

[TABLE]

where $\kappa$ is the flux-mean opacity (per unit mass) and $c$ is the speed of light; the numerical value equates $\kappa$ with the electron-scattering opacity for a fully-ionized solar-abundance mix. The equatorial value, $g_{\rm e}^{\rm Edd}$, is plotted in Fig. 3; it is exceeded by the actual equatorial gravity only for $M\gtrsim 18\text{M}_{\odot}$.101010We are aware that this has implications for the topology of the equipotential surface that are neglected in the present work. However, for traditional von Zeipel gravity darkening, there is no effect (Howarth, 1997), and it is only at near-critical rotation that the ELR formalism departs significantly from this. If anything, radiation-pressure effects will render our lower limits to mass more secure.

Also shown in Fig. 3 is a mean gravity, $\overline{g}=\int{g^{\ell}\,\text{d}a}/\int{\text{d}a}$, where $g^{\ell}$ is the local gravity and the integration is over the projected area. If we suppose that the observed $\log_{10}(g)$ values fall between the models’ equatorial and polar values then $M\lesssim 40\text{M}_{\odot}$; or if we speculate that the observed gravity can be identified with the mean gravity, then $20\lesssim M/\text{M}_{\odot}\lesssim 30$.

Table 1 includes detailed results for models at 15, 25, and 50$\text{M}_{\odot}$, representing what we consider to be the very extreme range of plausible masses if $P_{\rm rot}=1.78$ d, together with a ‘best guess’ central value. These models are intended to illustrate the sensitivity, or otherwise, of various parameters to the assumed mass.111111In principle, the synthetic spectra generated as part of the modelling could be used to gravity-sensitive lines, such as the wings of H$\gamma$, to better constrain the mass. Unfortunately, it isn’t possible to achieve convergence for line-blanketed, hydrostatic models for gravities less than $\log_{10}(g)\lesssim 3.5$ at the relevant temperatures; in our exoBush modelling, where necessary we used the lowest available gravity. The inferred inclination and equatorial velocity are particularly robust, with mean values averaged over all models, 15–50$\text{M}_{\odot}$, being $i=33{\aas@@fstack{\circ}}2\pm 1{\aas@@fstack{\circ}}8$, $v_{\rm e}=390\pm 16$ km s*-1*.

The implied axial inclination if $P_{\rm rot}=1.78$ d is unremarkable, but the equatorial velocity would be unprecedented for an early O I star; in the Galaxy, only a handful of near-main-sequence, late-O stars have comparable (projected) rotation velocities.121212 HDs 93521, 149757, and the ONn stars (Howarth & Smith, 2001; Walborn et al., 2011; Martins et al., 2015b).

4.5 Evolutionary status

The Hipparcos distance constrains the radius, hence luminosity, reasonably tightly, allowing us to locate $\zeta$ Pup quite precisely in the H–R diagram (Fig. 4), independently of its supergiant spectral classification. It is evident, simply by inspection, that our empirically inferred preferred mass range, $\sim$20–30$\text{M}_{\odot}$, is inconsistent with the evolutionary mass implied by the Brott et al. (2011) models, $M_{\rm evol}\sim 40$–50$\text{M}_{\odot}$ (rapidly-rotating$\rightarrow$non-rotating progenitor).

We can elaborate this inference with a Bonnsai analysis (Schneider et al. 2014131313The BONNSAI web service is available at

http://www.astro.uni-bonn.de/stars/bonnsai), built on the same evolutionary models. For a minimal set of observational constraints ($T_{\rm eff}=40\pm 1$ kK, $\log_{10}\left({L/\text{L}_{\odot}}\right)=5.60\pm 0.05$, $v_{\rm e}\sin{i}=213\pm 7$ km s*-1*, and a default set of priors), the single-star evolutionary tracks imply $M_{\rm evol}/\text{M}_{\odot}=42.4^{+3.4}_{-3.6}$ (Table 2, ‘Model 1’).

This minimal model has a predicted surface-helium abundance that is essentially solar, in conflict with observations; the sources listed in Section 4.1.2 give helium abundances by number in the range $y=0.14$–0.20, averaging $0.17\pm 0.02$ (s.d.; i.e., mass fraction $Y\simeq 0.41\pm 0.05$). A solution is still possible after imposing this additional constraint (‘Model 2’), with $M_{\rm evol}/\text{M}_{\odot}=37.0^{+3.8}_{-2.9}$.

However, imposing the further constraint of $v_{\rm e}=390\pm 16$ km s*-1* allows for no acceptable Bonnsai solutions (at any $Y$).

The discrepancies between empirical and single-star evolutionary masses support the proposal by van Rensbergen et al. (1996) that the runaway and rotational properties of $\zeta$ Pup are most readily understood in the context of previous binary interaction, and not single-star evolution. Although recent versions of this ‘Brussels scenario’ target high-mass solutions for $\zeta$ Pup, with M\gtrsim 60$$\text{M}_{\odot}(e.g., Vanbeveren 2012, Pauldrach et al. 2012; Ramiaramanantsoa et al. 2018), alternative channels can lead to undermassive (or overluminous) runaways (Vanbeveren & De Loore, 1994).

5 Starspots?

Non-radial pulsations and corotating starspots are the two most obvious candidates for the processes underpinning the 1.78-d photometric variability observed in $\zeta$ Pup (cf. $\S$3.4.1 of Ramiaramanantsoa et al. 2018 for a thorough discussion). Howarth & Stevens (2014) suggested that the physical origin of the signal may be pulsation associated with low-$\ell$ oscillatory convection modes, noting rough consistency with theoretical models by Saio (2011). Our revised physical parameters render the comparison problematic; Saio’s stability analysis was based on structures computed for standard single-star evolutionary tracks, which don’t explore the parameter space that now appears pertinent to $\zeta$ Pup ($T_{\rm eff}\simeq 40$ kK, $M\simeq 25\text{M}_{\odot}$, with potentially strong mixing in the scenario sketched by Vanbeveren & De Loore 1994).

By contrast, Ramiaramanantsoa et al. (2018) advocated a rotationally-modulated ‘hotspot’ interpretation of the 1.78-d period, and on that basis used light-curve inversion techniques to map the required surface-brightness distribution. Our updated understanding of the stellar geometry allows us to re-examine this question, and in particular to assess if line-profile variability offers a potential test of the rotational-modulation hypothesis.

5.1 Photometry

For an initial exploration, we took a representative model with $P_{\rm rot}=1.78$ d and $M=25\text{M}_{\odot}$ (cf. Table 1). We modified exoBush to include starspots that subtend a constant angular radius $\alpha_{\rm s}$ at the centre of mass (i.e., are approximately circular on the surface), at constant temperature $T_{\rm s}$. For specificity, we compare results with Ramiaramanantsoa et al.’s ‘Part IV’ photometric dataset, which has a moderately, but not exceptionally, large amplitude ($\sim$1%; Fig. 5) and their corresponding $i=33^{\circ}$ surface-brightness reconstruction (cf. panel 2 of their Fig. 10). This reconstruction has a rather simple hotspot geometry; we approximate their results with two $\sim$equatorial spots (colatitudes $\theta_{\rm s}=85^{\circ}$), separated by $158^{\circ}$ in longitude. While this spot model is not intended as a detailed ‘best fit’, it captures the essential characteristics of the observations at this epoch.

Based on the Ramiaramanantsoa et al. (2018) reconstruction, values of $\alpha_{\rm S}=10^{\circ}$, $T_{\rm s}=42.5$ kK were first employed; the results of the ($V$-band) photometric predictions are confronted with observations in Fig. 5. They substantially underestimate the observed amplitude; we believe this to be a straightforward consequence of the fact that the inversion approach is essentially a mathematical methodology intended simply to recover a surface-brightness distribution, while ours is a direct physical model.

To recover the amplitude observed in the BRITE-Constellation photometry requires spots that are either significantly larger or significantly hotter than initially assumed. We find that $\alpha_{\rm s}\simeq 26^{\circ}$ for $T_{\rm s}=42.5$ kK; $T_{\rm s}=50$ kK requires $\alpha_{\rm s}\simeq 14^{\circ}$.

Although these models reproduce the general characteristics of the photometry, we recall that it is always possible to construct a spot model capable of reproducing periodic, low-amplitude photometric variability; and a successful model fit (or inversion) based on the assumption of surface hotspots is not a proof of their existence, but is only a plausibility check. To emphasize this point, we have constructed a simple illustrative model with two cool spots, $T_{\rm s}=30$ kK, $\alpha_{\rm s}(1,2)=(17,20)^{\circ}$, separated by 180∘ in longitude. This ad hoc model matches the data at least as well as the hotspot calculations (Fig. 5).

5.2 Line-profile variability

A potentially more stringent test of the nature (and existence) of any starspots is offered by the spectroscopic line-profile variability that they should generate. An exploratory ‘proof of concept’ calculation is shown in Fig. 6, based on the $2\times(T_{\rm s}=50$ kK, $\alpha_{\rm s}=14^{\circ})$ spot model.

Cool spots generally give rise to (pseudo-)emission ‘bumps’ in spectrally resolved line profiles, essentially because they remove less flux at the projected spot velocity than in the continuum (e.g., Vogt & Penrod 1983). Bright spots will give rise to ‘absorption dips’ by the same reasoning, but only as long as the line equivalent widths do not change by a significant amount between the spot and adjacent unperturbed photosphere.

That is not the case for the model presented here. The assumed temperature contrast, $T_{\rm e}$ : $T_{\rm s}\sim 37$:50, is strong enough to result in substantial changes in line strength. For example, C iii $\lambda$5826 is a weak emission line in the intrinsic equatorial spectrum, but disappears entirely in the model’s high-temperature spot; similarly, He i $\lambda$5876 absorption weakens greatly. Moreover, the spots are hot, at low gravity, and viewed relatively far from normal incidence. In the hydrostatic, line-blanketed models used here, these circumstances are often accompanied by substantial decreases in the strengths of absorption lines, which may even go into emission through nLTE effects, as occurs for C iv $\lambda\lambda$5801, 5812. Consequently, when rectified and differenced with the unperturbed spectrum, the signature variability of hot spots in the models can be in the form of either absorption or emission features (Fig. 6).

Whether or not these modelled line-strength characteristics are quantitatively reliable, there is considerable diagnostic potential in the straightforward dynamical content of velocity- and temporally-resolved spectroscopy. Most importantly, Fig. 6 illustrates a specific discriminant between traditional, $i\simeq 90^{\circ}$, models and the low-inclination, $P_{\rm rot}\ =P_{\rm phot}$ alternative. For $\sim$equator-on configurations, longitudes more than $90^{\circ}$ from the central meridian are never visible, but at lower inclinations (as required by rotational modulation), features that occur in the hemisphere nearer the observer (i.e., ‘north’ of the equator) are visible beyond $\pm$0.25 in phase either side of transit, at submaximal velocity excursions, generating ‘

$\sim$

’-shaped features with red-to-blue ‘tails’ in the dynamic spectrum.

This is true even for the near-equatorial ($\theta_{\rm s}=85^{\circ}$), slightly extended spots used in the exploratory model. Spots with smaller colatitudes would be visible for greater fractions of the rotation period (giving relatively bigger tails to the

$\sim$

 tracks), and would also leave their signature over smaller ranges in velocity ($\Delta{v_{\rm s}}/v_{\rm e}\sin{i}\sim\sin\theta_{\rm s}$).

The predicted amplitudes of the spectroscopic spot signatures are on the order of $\sim$1%, which should be comfortably observable. Ramiaramanantsoa et al. (2018; their $\S$4.1.3) report recovering the 1.78-d signal in photospheric absorption lines, including C iv $\lambda\lambda$5801, 5812, but with “no obvious pattern”, which may be seen as weak evidence against rotational modulation. Other published time series also show only marginal evidence, at best, for any $\sim$1.8-d periodicity (Reid & Howarth, 1996; Berghöfer et al., 1996), although we have no way of knowing if the photometric signal was present at the time of those observations (cf. $\S$3).

6 Photospheric and stellar-wind variability

6.1 Corotating Interaction Regions and Discrete Absorption Components

Ultraviolet spectroscopy with the International Ultraviolet Explorer (IUE; Boggess et al. 1978) showed that ‘discrete absorption components’, or DACs, are a ubiquitous characteristic of early-type stars with strong winds. They are characterized by red-to-blue migration of features through the absorption components of P-Cygni profiles, with accelerations that are significantly slower than expected for the ambient outflow (e.g., Prinja & Howarth 1986, 1988; Kaper et al. 1996, 1999).

The slow acceleration, in particular, has prompted an interpretation in terms of corotating interaction regions (CIRs; Mullan 1984, 1986; Cranmer & Owocki 1996). This interpretation is bolstered by the (rather loose) anticorrelation between DAC recurrence timescales and $v_{\rm e}\sin{i}$ values for individual stars (Prinja 1988; Kaper et al. 1999; Howarth 2007). However, in most cases the DACs are not strictly periodic,141414An apparently periodic feature has been observed in UV P-Cygni profiles of the rapidly rotating B0.5 Ib star HD 64760, but this seems to be distinct from classical DACs (Prinja et al., 1995; Fullerton et al., 1997). so it remains open to debate as to whether the DACs are initiated directly in corotating (or very nearly corotating) photospheric features, or develop as an intrinsic property of radiatively-driven winds (cf., e.g., Martins et al. 2015a), notwithstanding that DAC-like behaviour has been traced to rather low velocities (Massa & Prinja, 2015).

6.2 CIRs, DACs, and spots

The phenomenological two-dimensional radiation-hydrodynamical CIR simulations by Cranmer & Owocki (1996; see also David-Uraz et al. 2017) employed the heuristic mechanism of photospheric bright spots to drive locally enhanced mass outflows, which they showed can generate DAC-like features (through velocity plateaux, rather than directly through density enhancements). Such spots may arise through subphotospheric convection zones (Cantiello & Braithwaite, 2011), and subsequent observational efforts have sought to identify corresponding spot-like surface features, and to associate them with DAC drivers (e.g., Ramiaramanantsoa et al. 2014).

$\zeta$ Pup is particularly well suited to such investigations, as it has the longest intensive UV-spectroscopy time series of any O star: a 16-d sequence of IUE observations in 1995 January (Howarth et al., 1995). These reveal a DAC recurrence timescale $T_{\rm DAC}=19.23{\pm}0.45$ hr at that epoch. Essentially the same timescale was found by Reid & Howarth (1996) in H$\alpha$, over the velocity range $\sim{-300}:-800$km s*-1*, in observations obtained in 1990.

Unfortunately for the ‘corotating spots’ hypothesis, the 1.78-d period is not commensurate with this DAC timescale; in particular, the duration of the IUE time series (hence precision of the timescale determination) rules out the possibility that $P_{\rm phot}=2\times T_{\rm DAC}$ with 4-$\sigma$ confidence.

Ramiaramanantsoa et al. (2018) sought to reconcile the superficially inconsistent timescales of the 1.78-d photometric signal and (twice) the $\sim$19-hr DAC recurrence period by speculating that $\zeta$ Pup may rotate differentially. In that case, a direct causal link between the mooted photospheric hot-spots and CIR formation could be salvaged if the spots occurred at higher, faster-rotating latitudes at the time of the IUE observations than at the epochs of the space-based photometry.151515The required differential rotation is antisolar – i.e., in the opposite sense to that observed in the Sun and solar-type stars, where the rotation period is shortest at the equator (e.g., Carrington 1860; Benomar et al. 2018). This speculation is bolstered by growing evidence for differential rotation in some stars bluewards of the granulation boundary (Balona & Abedigamba, 2016).

6.3 Differential rotation?

Direct testing of this speculation would require contemporaneous, extensive time series of both UV spectroscopy and photometry, which are unlikely to be available in the immediate future. However, differential rotation also has observable consequences for photospheric line profiles. To explore this, we performed calculations using a simple prescription for latitudinally differential angular rotation $\omega$ at colatitude $\theta$,

[TABLE]

(which in the spherical limit simplifies to the form commonly adopted for late-type stars, ${\omega(\theta)}/{\omega_{\rm e}}=1-k\cos^{2}\theta$). To estimate the $k$ parameter, we suppose that at the time of the IUE observations any spot features transited centrally (i.e., $\theta_{\rm s}\simeq i$), as suggested by the large covering factor of the DACs (i.e., the large fractional coverage of the projected stellar disk), and that their rotation period was $2T_{\rm DAC}$ (with an equatorial rotation period of $P_{\rm phot}$= 1.78 d); these assumptions lead to $k=-0.16$.

Three sets of model profiles are shown in Fig. 7, based on general physical parameters from columns 3 and 5 of Table 1. We stress that these are ad hoc, ab initio model calculations, and are not, in any sense, fits to observations (which would properly entail exploration of a wider parameter space). Nevertheless, the comparison with observed profiles is of some interest; the traditional, $i\simeq 90^{\circ}$ ($P_{\rm rot}=3.35$ d) model provides a reasonably satisfactory match to the C iv absorption profiles, while the $P_{\rm rot}=1.78$ d model fares rather less well. The discrepancies with the differentially rotating model are large enough to cast doubt on the underpinning speculation, and hence on the proposal that DACs in $\zeta$ Pup are directly driven by the same phenomenon that is responsible for the photometric signal.

We mention two further practical difficulties that challenge a model whereby DACs are the result of CIRs driven directly by corotating photospheric hot spots:

   (i) It appears likely that all early-type stars with strong winds exhibit DACs (e.g., Howarth & Prinja 1989) and yet there is currently little evidence for ubiquitous periodic photospheric spot activity (or for strictly periodic DAC activity). In the specific case of $\zeta$ Pup, there is a clear record of ‘normal’ DAC behaviour in 1989 April (Prinja et al., 1992), while contemporaneous photometry gives no indication of significant variability on a commensurate period (Balona 1992; $\S$3).

   (ii) The DAC covering factors are large ($\gtrsim$0.5; e.g., Howarth & Smith 1995; Massa & Prinja 2015). If we suppose that DACs were present in $\zeta$ Pup at the epochs of the 21st-century space photometry (which, though now untestable, seems likely, given their $\sim$universal occurence), it is not obvious how this can be reconciled with an origin in CIRs originating from relatively small, $\sim$equatorial spots viewed at low axial inclination.

6.4 Other timescales

Part of the motivation that led Ramiaramanantsoa et al. (2018) to seek an association of $P_{\rm phot}$ with $T_{\rm DAC}$ was their discovery that the 1.78-d period is not associated exclusively with the photosphere, but can additionally be traced in the He ii $\lambda$4686 emission line at velocities out to $\sim{\pm}400$ km s*-1*. It is apparent, therefore, that there is some connection between activity in the photosphere and the base of the stellar wind.

A similar association was reported by Reid & Howarth (1996; data obtained 1992), who found an 8.5-hour signal in photospheric lines, with blue-to-red propagation.161616This period was first identified by Baade (1986), in data obtained in 1984. He proposed an interpretation in terms of sectoral-mode non-radial pulsations ($\ell=|m|=2$; Baade 1988, Reid & Howarth 1996). The period has not been recovered in other datasets of comparable, or better, quality (cf. Howarth & Stevens 2014; Ramiaramanantsoa et al. 2018). At the risk of invoking arbitrary numerology, we note the coincidence that $P_{\rm phot}=1.78$ d is exactly 5$\times$ the 8.54-hr period reported by Reid & Howarth. They identified the same period in H$\alpha$ emission at velocities more negative than $-280$ km s*-1*, moving red to blue – further evidence for a stellar-wind signature of photospheric activity.

However, in that case the 8.5-hr signal co-existed with a 19.2-hr signal (matching $T_{\rm DAC}$), the latter again featuring red-to-blue migration, detected at outflow velocities from $\sim$300 to 800 km s*-1*. It therefore appears possible that the low-velocity wind may respond to photospheric drivers (of whatever nature) without that response necessarily propagating directly to the DACs in the high-velocity wind.

7 Summary & conclusion

We have argued that the Hipparcos astrometry for $\zeta$ Pup is reliable, and hence that $d=332\pm 11$ pc ($\S$2.1). With this distance, we modelled the basic physical characteristics under two extreme assumptions:

(i) $i\simeq 90^{\circ}$ ($P_{\rm rot}\simeq 3.3$ d, $v_{\rm e}\simeq 213$ km s*-1*, $\omega_{\rm e}/\omega_{\rm crit}\simeq 0.6$), and

(ii) $P_{\rm rot}=P_{\rm phot}$ ($i\simeq 33^{\circ}$, $v_{\rm e}\simeq 390$ km s*-1*, $\omega_{\rm e}/\omega_{\rm crit}\simeq 0.9$), where $P_{\rm phot}$= 1.78 d is the period observed in 21st-century broad-band, space-based photometry.

In either case, $R_{*}\simeq 13\text{R}_{\odot}$, $M\simeq 25\text{M}_{\odot}$ (cf. Table 1 and Section 4.4 for details). These characteristics are not consistent with single-star evolutionary tracks. As proposed by van Rensbergen et al. (1996) and Vanbeveren (2012), binary (or multiple-star) interaction earlier in $\zeta$ Pup’s lifetime is probably required, and is then implicated in its runaway status, and in its exceptionally rapid rotation.

The 1.78-d photometric variability was of smaller amplitude (or absent) in Hipparcos photometry than more-recent datasets ($\S$3), confirming previous indications that secular changes occur.

We have reviewed the possibility that the 1.78-d signal arises from rotational modulation of photospheric starspots (i.e., that $P_{\rm phot}=P_{\rm rot}$). Available modelling of the photometric data is incapable of discriminating between hot and cool spots (or other mechanisms, such as non-radial pulsation). Time-series spectroscopy resolving the line profiles should afford a reasonably clean test of the corotating-spot hypothesis (more precisely, of the implied low axial inclination; $\S$5.2), though existing data do not show the expected signature. We consider that the origin of the 1.78-d period remains open to question (particularly given that a separate, 2.56-d signal was present in the Hipparcos photometry, and was similarly attributed to rotational modulation).

We have re-examined the putative association between $P_{\rm phot}$, $P_{\rm rot}$, CIRs, and DACs. A direct causal association between DACs and corotating interaction regions driven at $P_{\rm phot}$ would require an ad hoc mechanism of differential rotation coupled to latitudinally mobile hotspots. Exploratory line-profile calculations afford a reasonable match to observations for $i=90^{\circ}$ ($P_{\rm rot}\simeq 3.3$ d); lower-inclination, differentially rotating models give poorer agreement ($\S$6.3). We conclude that a compelling case for DACs being the result of CIRs driven directly by corotating photospheric hotspots has yet to be made for $\zeta$ Pup.

acknowledgements

We thank Tony Moffat, Yaël Nazé, Stan Owocki, Raman Prinja, Tahina Ramiaramanantsoa, Dominic Reeve, Dany Vanbeveren, and Gregg Wade for useful correspondence/conversations.