# BRITE-Constellation high-precision time-dependent photometry of the early O-type supergiant ζ Puppis unveils the photospheric drivers of its small- and large-scale wind structures

Tahina Ramiaramanantsoa, Anthony Moffat, Robert Harmon, Richard Ignace, Nicole St-Louis, Dany Vanbeveren, Tomer Shenar, Herbert Pablo, Richard Richardson, Ian Howarth, et al.

► To cite this version:

Tahina Ramiaramanantsoa, Anthony Moffat, Robert Harmon, Richard Ignace, Nicole St-Louis, et al.. BRITE-Constellation high-precision time-dependent photometry of the early O-type supergiant ζ Puppis unveils the photospheric drivers of its small- and large-scale wind structures. Monthly Notices of the Royal Astronomical Society, 2018, 473 (4), pp.5532-5569. 10.1093/mnras/stx2671. hal-02106861

HAL Id: hal-02106861

https://hal.science/hal-02106861v1

Submitted on 22 Nov 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# BRITE-Constellation high-precision time-dependent photometry of the early O-type supergiant ζ Puppis unveils the photospheric drivers of its small- and large-scale wind structures

Tahina Ramiaramanantsoa,¹,²★ Anthony F. J. Moffat,¹,² Robert Harmon,³

Richard Ignace,⁴ Nicole St-Louis,¹,² Dany Vanbeveren,⁵ Tomer Shenar,⁶

Herbert Pablo,¹,² Noel D. Richardson,⁷ Ian D. Howarth,⁸ Ian R. Stevens,⁹

Caroline Piaulet,¹ Lucas St-Jean,¹ Thomas Eversberg,¹⁰ Andrzej Pigulski,¹¹

Adam Popowicz,¹² Rainer Kuschnig,¹³ Elżbieta Zocłońska,¹⁴ Bram Buysschaert,¹⁵,¹⁶

Gerald Handler,¹⁴ Werner W. Weiss,¹³ Gregg A. Wade,¹⁷ Slavek M. Rucinski,¹⁸

Konstanze Zwintz,¹⁹ Paul Luckas,²⁰ Bernard Heathcote,²¹ Paulo Cacella,²²

Jonathan Powles,²³ Malcolm Locke,²⁴ Terry Bohlsen,²⁵ André-Nicolas Chené,²⁶

Brent Miszalski,²⁷,²⁸ Wayne L. Waldron,²⁹ Marissa M. Kotze,²⁷,²⁸ Enrico J. Kotze²⁷ and Torsten Böhm³⁰,³¹

*Affiliations are listed at the end of the paper*

Accepted 2017 October 11. Received 2017 September 22; in original form 2017 May 23

## ABSTRACT

From 5.5 months of dual-band optical photometric monitoring at the 1 mmag level, *BRITE-Constellation* has revealed two simultaneous types of variability in the O4I(n)fp star ζ Puppis: one single periodic non-sinusoidal component superimposed on a stochastic component. The monoperiodic component is the 1.78-d signal previously detected by *Coriolis/Solar Mass Ejection Imager*, but this time along with a prominent first harmonic. The shape of this signal changes over time, a behaviour that is incompatible with stellar oscillations but consistent with rotational modulation arising from evolving bright surface inhomogeneities. By means of a constrained non-linear light-curve inversion algorithm, we mapped the locations of the bright surface spots and traced their evolution. Our simultaneous ground-based multisite spectroscopic monitoring of the star unveiled cyclical modulation of its He II λ4686 wind emission line with the 1.78-d rotation period, showing signatures of corotating interaction regions that turn out to be driven by the bright photospheric spots observed by *BRITE*. Traces of wind clumps are also observed in the He II λ4686 line and are correlated with the amplitudes of the stochastic component of the light variations probed by *BRITE* at the photosphere, suggesting that the *BRITE* observations additionally unveiled the photospheric drivers of wind clumps in ζ Pup and that the clumping phenomenon starts at the very base of the wind. The origins of both the bright surface inhomogeneities and the stochastic light variations remain unknown, but a subsurface convective zone might play an important role in the generation of these two types of photospheric variability.

**Key words:** techniques: photometric – techniques: spectroscopic – stars: massive – stars: rotation – starspots – supergiants – stars: winds, outflows.

## 1 INTRODUCTION

Spaceborne high-precision photometric monitoring of stars as a means to probe their intrinsic light variations and trace back the

<footnote>* E-mail: tahina@astro.umontreal.ca</footnote>
underlying physical mechanisms has only emerged over the past two decades (Buzasi 2000; Buzasi et al. 2000). Particularly, this practice is even less frequent on hot luminous massive O-type stars (with $M_{\text{init}} \gtrsim 20 M_{\odot}$), partly because of their scarcity which is reflected in the decreasing initial stellar mass function (IMF) following a power law of index −2.35 for masses between 1.25 and 150 $M_{\odot}$ (Salpeter 1955; Bastian, Covey & Meyer 2010). Nevertheless, with their strong radiation-driven stellar winds and their termination as supernovae, O stars and their descendant Wolf–Rayet stars are important drivers of the chemical enrichment of galaxies and the Universe. Understanding the physical origin of their intrinsic variations may provide constraints not only on their internal properties (e.g. convective core overshoot, possible existence of a subsurface convection zone, and magnetic fields) but also their global wind properties, and the complex physics of their photosphere–wind interface.

## 1.1 Space photometry of O stars

So far a total of 18 O stars have been monitored through space photometry (table 3 of Buysschaert et al. 2015, adding the recently reported photometric observations of the δ Orionis system by Pablo et al. 2015, the ι Orionis binary system by Pablo et al. 2017, the magnetic hot supergiant ζ Orionis Aa by Buysschaert et al. 2017, and the late O-type supergiant HD 188209 by Aerts et al. 2017). The Wide-field InfraRed Explorer (WIRE, Buzasi 2000) was the pioneer of this practice by including the late O-type subgiant ζ Oph in its sample of 90 targets (Bruntt 2007). Then followed contributions from the Microvariability and Oscillations of Stars microsatellite (MOST, Walker et al. 2003), the Convection Rotation and planetary Transits satellite (CoRoT; Baglin et al. 2006; Auvergne et al. 2009), the Solar Mass Ejection Imager (SMEI, Eyles et al. 2003) instrument onboard the Coriolis satellite, and the Kepler spacecraft (Borucki et al. 2010; Howell et al. 2014). Various types of intrinsic variations were reported from these space-based photometric observations, namely β Cep-like pulsations for the late O-type stars (ζ Oph, Walker et al. 2005, Howarth et al. 2014; HD 46202: Briquet et al. 2011; HD 256035: Buysschaert et al. 2015), one case of the presence of p modes with finite lifetimes analogous to solar-like oscillations (HD 46149, Degroote et al. 2010), two cases of tidally influenced pulsations (δ Ori Aa, Pablo et al. 2015; ι Ori, Pablo et al. 2017), one case of a possible oscillatory convection g mode (ζ Pup, Howarth & Stevens 2014), six cases of possible internal gravity waves (IGWs: HD46223, HD46150, HD 46966, Blomme et al. 2011; Aerts & Rogers 2015; HD44597, HD 255055, Buysschaert et al. 2015; HD 188209, Aerts et al. 2017) and five cases of possible rotational modulation (HD 46149, Degroote et al. 2010; HD 46223, Blomme et al. 2011; ξ Per, Ramiaramanantsoa et al. 2014; EPIC 202060097, EPIC 202060098, Buysschaert et al. 2015). Notably, some stars show simultaneously different types of intrinsic variations. With its sample of 36 O stars brighter than $V = 6$, the recently launched, commissioned and now operational nanosatellites forming BRight Target Explorer (BRITE-Constellation, Weiss et al. 2014; Pablo et al. 2016) will almost triple the current numbers and improve the view on this broad range of O star intrinsic variability.

## 1.2 ζ Puppis

Standing out from the aforementioned restricted sample of O-type stars is the hot, bright, rapidly rotating, single, runaway supergiant

<table><caption>**Table 1.** Stellar parameters for ζ Puppis (HD 66811). The parameter β denotes the exponent of the usual β-law velocity generally adopted in hot stellar winds, while $f_{\infty}$ denotes the asymptotic value of the clumping filling factor at the outer boundary of the stellar wind.</caption><thead><tr><th>Parameter</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>Spectral type</td><td>O4I(n)fp</td><td>(1)</td></tr><tr><td>V</td><td>2.256 ± 0.019</td><td>(2)</td></tr><tr><td>B − V</td><td>−0.276 ± 0.012</td><td>(2)</td></tr><tr><td>U − B</td><td>−1.108 ± 0.013</td><td>(2)</td></tr><tr><td>v<sub>e</sub> sin i</td><td>(km s<sup>−1</sup>)</td><td>(3)</td></tr><tr><td>v<sub>rad</sub></td><td>(km s<sup>−1</sup>)</td><td>(4)</td></tr><tr><td>v<sub>tan</sub></td><td>(km s<sup>−1</sup>)</td><td>(5)</td></tr><tr><td>Distance</td><td>(kpc)</td><td>(6)</td></tr><tr><td>Age</td><td>(Myr)</td><td>(6)</td></tr><tr><td>log (L/L<sub>☉</sub>)</td><td>(cgs dex)</td><td>(6)</td></tr><tr><td>T<sub>eff</sub></td><td>(kK)</td><td>(6)</td></tr><tr><td>log g</td><td>(cm s<sup>−2</sup>)</td><td>(6)</td></tr><tr><td>R</td><td>(R<sub>☉</sub>)</td><td>(7)</td></tr><tr><td>M</td><td>(M<sub>☉</sub>)</td><td>(6)</td></tr><tr><td>M&#775;</td><td>(M<sub>☉</sub> yr<sup>−1</sup>)</td><td>(6)</td></tr><tr><td>v<sub>∞</sub></td><td>(km s<sup>−1</sup>)</td><td>(6)</td></tr><tr><td>β</td><td></td><td>(6)</td></tr><tr><td>f<sub>∞</sub></td><td>0.05</td><td>(6)</td></tr></tbody></table>

ζ Puppis (Table 1). Owing to its brightness, proximity and status as a single star, ζ Pup has become a key object for understanding the properties of O-type stars and hot stellar winds. Therefore, it is not surprising why several observational studies have been conducted on ζ Pup across a wide range of wavelength: in the X-ray (e.g. Hervé, Rauw & Nazé 2013; Nazé, Oskinova & Gosset 2013, and references therein), in the ultraviolet (UV, Prinja et al. 1992; Howarth, Prinja & Massa 1995), in the optical (Davis et al. 1970; Reid & Howarth 1996; Howarth & Stevens 2014), as well as in the radio and submillimetre (Blomme et al. 2003). Amongst these studies, time-resolved spectroscopic and photometric observations have unveiled a number of variabilities on time-scales ranging from hours to a couple of days (table 1 in Howarth & Stevens 2014).

The first variability reported on ζ Pup came from the monitoring of the central absorption reversal of its Hα wind emission-line profile which was seen to vary non-sinusoidally with a period of 5.075 ± 0.003 d (Moffat & Michaud 1981). Subsequent ground-based optical photometric observations in the Strömgren b filter (Balona 1992), as well as UV spectroscopic studies (Howarth et al. 1995) revealed periods compatible with this one (~5.2 d). Up to this time, this 5.1 d period was attributed to the rotation period of the star, which, combined with the measured projected rotational velocity and the radius (Table 1), implies that the star is seen almost equator on.

Also, photospheric line-profile variations (LPVs) that were best interpreted as non-radial pulsations (NRPs) with $l = -m = 2$ were detected in ζ Pup, occurring as bumps moving redward across some optical He I, He II, N IV and C IV lines of the star with a period of 8.54 ± 0.05 h (Baade 1986; Reid & Howarth 1996). Additionally, a possibility for the presence of low-order oscillatory convection modes with a period of 1.780938 ± 0.000093 d and semi-amplitude of 6.9 ± 0.3 mmag has been reported from four contiguous years (2003–2006) of Coriolis/SMEI photometric monitoring of the star (Howarth & Stevens 2014).

<footnote>(1) Galactic O-Star Catalog: Sota et al. (2014), (2) Maíz-Apellániz et al. (2004), (3) Howarth et al. (1997), (4) Gontcharov (2006), (5) Moffat et al. (1998), (6) Bouret et al. (2012), (7) Martins et al. (2005).</footnote>
From the point of view of wind variability, $\zeta$ Pup is often used as a laboratory to study the time-dependent properties of wind structures in O-type stars. Such wind structures are classified in two distinct categories with respect to their spatial extent: small-scale density inhomogeneities dubbed clumps; and large-scale spiral-like structures corotating with the stellar photosphere, the so-called corotating interaction regions (CIRs). Both of these phenomena are now believed to be essentially ubiquitous in O stars (Puls, Vink & Najarro 2008).

Observational evidence for the presence of clumps in hot luminous star winds was initially detected in stars in the He-burning, He-rich Wolf–Rayet phase (Schumann & Seggewiss 1975; Moffat et al. 1988), a late stage in the evolution of O-type stars (e.g. Lamers et al. 1991). In the case of O stars, a phenomenological model consisting of a population of radiatively driven blobs in the wind of $\zeta$ Pup was first proposed to explain its observed X-ray luminosity (Lucy & White 1980). Later $\zeta$ Pup became the first O-type star confirmed to have a clumpy wind: spectroscopic monitoring of its wind-formed He II $\lambda$4686 line profile revealed excess bumps at the ~5 per cent level of the line intensity, propagating from the centre towards the edges of the line, which is a signature of the kinematics of clumps randomly appearing in the inner stellar wind (Eversberg, Lépine & Moffat 1998). Subsequent studies confirmed these interpretations and detected similar behaviour in other hot stars (Lépine & Moffat 2008).

The physical origin of these clumps still remains to be determined, with two coexistent paradigms to be considered: the first scenario involves line-driving instability that could generate inverse shocks locally compressing wind material to form clumps in a less dense interclump medium (Lucy & Solomon 1970; Owocki, Castor & Rybicki 1988; Feldmeier 1995; Dessart & Owocki 2003), whereas the second one involves a contribution from IGWs that could be excited within a subsurface convection zone and induce the formation of clumps near the photosphere–wind interface (Cantiello et al. 2009) which could then be subsequently enhanced by the first mechanism further out in the stellar wind.

Besides that, the absorption troughs of unsaturated P Cygni profiles of UV resonance lines of O-type stars are the sites for the appearance and blueward propagation of spectral features called discrete absorption components (DACs), a more general form of features initially called narrow absorption components (Prinja 1988; Howarth & Prinja 1989; Howarth 1992; Howarth et al. 1995; Kaper et al. 1999). These features are best interpreted as the spectroscopic manifestations of the presence of CIRs in the winds of O stars (Mullan 1984), such as those that are also found in the solar wind (Hundhausen 1972). $\zeta$ Pup was one of the three program stars observed during extensive UV spectroscopic monitoring with the *International Ultraviolet Explorer* (IUE) (Massa et al. 1995; Howarth et al. 1995). A mean DAC recurrence period of 19.23 ± 0.45 h was measured from the analysis of their propagation in the Si IV $\lambda\lambda$1394, 1403 doublet line of the star. With a 5.1 d rotation period as was believed at that time, $\zeta$ Pup would exhibit on average 5–6 DACs per rotation cycle (fig. 1 in Howarth et al. 1995), which is much more than the average of two dominating DACs per rotation cycle observed in other O stars (Kaper et al. 1999).

As for the origin of these CIRs, the canonical hydrodynamical model in O stars involves photospheric perturbations such as corotating bright spots or NRPs which initially enhance the local stellar wind speed, generating wind plasma of a different speed that interacts with the ambient wind to produce compressed material shaped as large spiral arms as they corotate with the star (Cranmer

& Owocki 1996). Velocity plateaus forming ahead of the CIR arms are seen as propagating DACs in the absorption troughs of UV P Cygni resonance lines such as C IV $\lambda$1548 or Si IV $\lambda\lambda$1394, 1403. Observational evidence for the presence of bright spots on an O star with possible link to its known DACs recurrence period recently emerged from MOST photometric monitoring of the single runaway mid-O-type giant $\xi$ Per (Ramiaramanantsoa et al. 2014). This discovery raised a motivation to determine if the presence of bright photospheric spots is universal among O stars (like the NAC/DAC phenomenon) and how they are linked to CIR/DAC activities in the stellar wind. Primarily driven by this purpose, we conducted a coordinated optical observing campaign on $\zeta$ Pup spanning ~5.5 months, consisting of space-based photometric monitoring with BRITE-Constellation during an observing run on the Vela/Puppis field (2014 December 11–2015 June 02), and contemporaneous multisite ground-based spectroscopic monitoring of the He II $\lambda$4686 line profile of the star. Our BRITE observations probe light variations at the level of the stellar photosphere in two passbands, not achievable by any previous space-based mission, while our spectroscopic observations of the He II $\lambda$4686 wind emission line allow us to trace any possible activities related to CIRs and clumps in the inner stellar wind.

## 2 OBSERVATIONS

### 2.1 High-precision space photometry

#### 2.1.1 BRITE-Constellation photometry

BRITE-Constellation (Weiss et al. 2014; Pablo et al. 2016) as originally planned consists of a network of six nanosatellites each housing a 35-mm format KAI-11002M CCD photometer fed by a 30-mm diameter f/2.3 telescope through either a blue filter (390–460 nm) or a red filter (545–695 nm): BRITE-Austria (BAb), Uni-BRITE (UBr), BRITE-Heweliusz (BHr), BRITE-Lem (BLb), BRITE-Toronto (BTr) and BRITE-Montréal (BMb), the last letter of the abbreviations denoting the filter type ('b' for blue and 'r' for red). With the exception of BMb which failed to detach from the upper stage module of its Dnepr launch vehicle, all the satellites were launched into Low Earth Orbits of orbital period of the order of 100 min, commissioned and now fully operational. With a ~24° × 20° effective field of view, each component of BRITE-Constellation performs simultaneous monitoring of 15–30 stars brighter than $V \simeq 6$. A given field is observed typically over a ~6 month time base and, as far as possible, at least two satellites equipped with different filters are set to monitor the field to ensure dual-band observations.

The seventh field monitored by BRITE was the Vela/Puppis field, for which photometry of 32 stars was extracted, $\zeta$ Pup being the sole O-type star (not counting the primary component of the WC+O system $\gamma^2$ Velorum, as the binary is not resolved by BRITE, the BRITE detector pixel size being ~27.3 arcmin). Observations of this field were performed by BAb, BTr and BHr (Table 2) during ~5.5 months between 2014 December 11 and 2015 June 02 (HJD 2, 457, 003.40 – 176.31). Short 1 s exposures were taken at a median cadence of 15.3 s during ~1 per cent–30 per cent of each ~100 min BRITE orbit, the remaining time unused due to stray light interferences, blocking by the Earth and limited data download capacity. Observations were performed in stare mode for all satellites and then switched to chopping mode (Pablo et al. 2016; Popowicz 2016; Popowicz et al. 2017) for BAb and BTr for the second half of the observing run. Onboard stacking of three consecutive frames was
<table><caption>**Table 2.** BRITE photometry of ζ Pup – journal of observations. $N_s$ indicates the number of frames stacked onboard. Quantities listed in columns 6–9 were assessed at post-correlation stage: $N_{p, tot}$ is the total number of data points, $N_{p, orb}$ is the median number of points per orbit (along with the range of values), $t_{cont}$ is the median contiguous time per orbit during which observations were performed (along with the range of values) and $\sigma_{rms}$ is the quadratic mean value of mean standard deviations per orbit assessed from linear fits of the flux values within the orbits.</caption><thead><tr><th>Satellite</th><th>Start date<br>(HJD-2451545)</th><th>End date</th><th>Observing<br>mode</th><th>$N_s$</th><th>$N_{p, tot}$</th><th>$N_{p, orb}$</th><th>$t_{cont}$<br>(min)</th><th>$\sigma_{rms}$<br>(mmag)</th></tr></thead><tbody><tr><td rowspan="2">BRITE-Austria</td><td>5458.402</td><td>5527.991</td><td>Stare</td><td>1</td><td>16198</td><td>25[4–48]</td><td>6.9[0.8–12.4]</td><td>1.41</td></tr><tr><td>5527.920</td><td>5621.152</td><td>Chopping</td><td>1</td><td>22981</td><td>30[6–40]</td><td>11.6[2.3–13.6]</td><td>1.19</td></tr><tr><td rowspan="2">BRITE-Toronto</td><td>5465.758</td><td>5535.193</td><td>Stare</td><td>3</td><td>15744</td><td>23[4–30]</td><td>16.8[2.2–22.0]</td><td>1.36</td></tr><tr><td>5535.135</td><td>5624.967</td><td>Chopping</td><td>1</td><td>54135</td><td>56[8–60]</td><td>15.1[2.3–15.6]</td><td>0.75</td></tr><tr><td>BRITE-Hewelusz</td><td>5553.059</td><td>5631.313</td><td>Stare</td><td>1</td><td>44747</td><td>78[8–161]</td><td>15.9[2.1–30.1]</td><td>1.40</td></tr></tbody></table>

performed only for the first three setups of the *BTr* observations. Raw light curves were extracted using the reduction pipeline for *BRITE* data which also includes corrections for intrapixel sensitivity (Popowicz et al. 2017). Then post-reduction decorrelations with respect to instrumental effects due to CCD temperature variations, centroid position and orbital phase were performed on each observational setup for each satellite according to the method described by Pigulski et al. (2016), and flux variations due to changes in point spread function shape as a function of temperature were performed according to the method described by Buysschaert et al. (2017). In the resulting final decorrelated light curves, we do not notice any obvious variations on time-scales shorter than the *BRITE* orbital period ($t \lesssim 100$ min) that could be qualified as intrinsic to the star rather than pure instrumental noise. Therefore, in order to gain in precision we calculated satellite-orbital mean fluxes to create the final light curves in the two filter bands. We note that we also perform removal of outliers during the decorrelation process, such that it is reasonable to adopt orbital mean fluxes instead of median fluxes or trimmed mean fluxes or mean flux values within 1/4 and 3/4 quartiles. Then, to extract the orbital means, we tested two different methods: (1) a simple average of the flux values taken within an orbit, and (2) an average taken to be the mid-point of a linear fit of the flux values within an orbit. We noticed no significant difference between the resulting root mean square (rms) values of the mean standard deviations obtained from the two methods, the second one being only slightly better. Therefore, we adopted the second method to generate the final binned light curves in the two filter bands which we use to extract information on the intrinsic variability of the star (Fig. 1).

### 2.1.2 Coriolis/SMEI photometry

We also use the archival light curves of ζ Pup recorded by Coriolis/SMEI during its 2003–2006 seasonal observing runs, published by Howarth & Stevens (2014). The characteristics of the Coriolis/SMEI light curves of ζ Pup are described in detail in section 2.1 of Howarth & Stevens (2014).

## 2.2 Ground-based multisite optical spectroscopy

In parallel with the *BRITE* run, we conducted multisite ground-based spectroscopy of ζ Pup focusing on the He II λ4686 line in order to monitor possible signatures of CIR/DAC activities in the inner stellar wind. The temporal coverage of the ground-based observations is depicted in Fig. 1 and the characteristics of the observations are listed in Table 3.

### 2.2.1 CTIO-SMARTS 1.5 m/CHIRON

We acquired a total of 415 spectra of ζ Pup spanning 4578–8762 Å with the CHIRON échelle spectrograph (Tokovinin et al. 2013) mounted on the Cerro Tololo Inter-American Observatory’s (CTIO) 1.5 m telescope operated by the Small and Moderate Aperture Research Telescope System (SMARTS) consortium (NOAO Proposal IDs: 2014B-0082, 2015A-0002; PI: T. Ramiaramantsoa). We adopted a cadence of 2–6 visits per night on ζ Pup (depending on the observability of the star), between 2014 November 28 and 2015 June 20. Each visit consists of 4–6 consecutive 15 s subexposures taken in slicer mode ($R \sim 80\,000$) that we stacked and rebinned down to $R \sim 24\,000$ so that one single spectrum has a $S/N \sim 725$ pixel$^{-1}$ in the continuum close to the He II λ4686 line. In order to correct for the instrumental blaze function, we normalized the extracted one-dimensional wavelength-calibrated spectra of ζ Pup through a division by cubic spline fits to each of the orders of the spectrum of the B9.5V star HR4468 ($V = 4.7$) obtained during the campaign with the same instrument configurations (slicer mode, fixed cross-disperser position). We chose HR4468 for this purpose as it is relatively close to ζ Pup, relatively bright enough and does not have many strong metal lines. Fig. 2 illustrates the normalized, unbinned ($R \sim 80\,000$) CTIO-SMARTS 1.5 m/CHIRON mean spectrum of ζ Pup over the entire campaign with identification of its strongest spectral lines.

### 2.2.2 SAAO 1.9 m/GIRAFFE

We collected 24 optical spectra of ζ Pup during 14 nights between 2015 January 28 and 2015 February 10 at the South-African Astronomical Observatory (SAAO) with the Grating Instrument for Radiation Analysis with a Fibre Fed Échelle (GIRAFFE) spectrograph hosted in the coudé chamber of the 1.9 m Grubb Parsons telescope (PI: T. Ramiaramantsoa). Given the two available cross-disperser prisms optimized for the blue domain (3791–5459 Å) and the red domain (4094–9397 Å) with typical $R \sim 39\,000$, we chose to use the red prism in order to better match with the CTIO-SMARTS 1.5 m/CHIRON wavelength coverage. Camera flat-fields, fibre flat-fields, bias frames and arc exposures were obtained in the usual way and the data reduced using the PYTHON-based pipeline for GIRAFFE data extraction, INDLULAMITHI. Exposure times range from 4 to 12 min depending on the weather conditions and the air mass of the star. We adopted a cadence of 2–5 visits per night well spread in time, alternating with visits on another primary target ($\gamma^2$ Velorum) and the B9.5V star HR4468 which was used to correct for the instrumental blaze function in the spectra of the two primary targets. Each visit on ζ Pup consists of 4–6 consecutive spectra that
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two-colour time series *BRITE* photometry of ζ Pup binned over each *BRITE* orbit. The observations obtained from *BAb*, *BTr* and *BHR* are in blue, red and dark cyan, respectively. For better visibility, an offset of 50 mmag was added to the light curve in the blue filter. The temporal coverage of the contemporaneous ground-based spectroscopy is also depicted as small diamonds below the *BRITE* light curves.</figcaption></figure>

we later stacked in order to get a S/N ~ 370 in the continuum near He II λ4686.

### 2.2.3 SASER

Optical spectra covering the He II λ4686 region were obtained by amateur observatories members of the Southern Astro Spectroscopy Email Ring (SASER) located in Australia (Shenton Park observatory, Domain Observatory, Latham Observatory, Mirranook Observatory), New-Zealand (R. F. Joyce Observatory) and Brazil

(Dogsheaven Observatory). Data reduction (bias, flat-field, dark and sky background correction) and extraction of one-dimensional wavelength-calibrated unnormalized spectra were performed with the Integrated Spectrographic Innovative Software (ISIS) and the MAXIM DL software. Then, by means of the Image Reduction and Analysis Facility (IRAF,<fnref n="1" /> Tody 1986, 1993) software, we performed

<footnote n="1">¹ IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.</footnote>
<table><caption>**Table 3.** Basic characteristics of the spectroscopic observations involved in our campaign. The facilities are ordered according to East longitude, starting from the closest to the Greenwich meridian. Amateur facilities members of the SASER (http://saser.wholemeal.co.nz) at the time of the campaign are indicated by a star. Estimates of the typical S/N per spectrum were taken in a portion of continuum in the closest vicinity of He II λ4686.</caption><thead><tr><th>Observatory</th><th>PI/Observer</th><th>East longitude</th><th>Telescope</th><th>Instrument</th><th>R</th><th>λ range (Å)</th><th>N<sub>sp</sub></th><th>S/N</th></tr></thead><tbody><tr><td>SAAO</td><td>T. Ramiaramanantsoa</td><td>20°48′37″</td><td>1.9 m</td><td>GIRAFFE</td><td>39000</td><td>4249–6513</td><td>24</td><td>370</td></tr><tr><td>Shenton Park *</td><td>P. Luckas</td><td>115°48′54″</td><td>0.35 m</td><td>Lhires III</td><td>8364</td><td>4579–4729</td><td>257</td><td>590</td></tr><tr><td>Domain*</td><td>B. Heathcote</td><td>144°59′21″</td><td>0.28 m</td><td>Lhires III</td><td>8529</td><td>4590–4738</td><td>106</td><td>700</td></tr><tr><td>Latham*</td><td>J. Powles</td><td>149°01′41″</td><td>0.25 m</td><td>Spectra L200</td><td>8774</td><td>4549–4879</td><td>65</td><td>500</td></tr><tr><td>Mirranook*</td><td>T. Bohlsen</td><td>151°30′35″</td><td>0.28 m</td><td>Spectra L200</td><td>6361</td><td>4484–5022</td><td>1</td><td>320</td></tr><tr><td>R. F. Joyce*</td><td>M. Locke</td><td>172°20′59″</td><td>0.4 m</td><td>Spectra L200</td><td>8432</td><td>4453–4809</td><td>22</td><td>725</td></tr><tr><td>CTIO</td><td>T. Ramiaramanantsoa</td><td>290°48′23″</td><td>1.5 m</td><td>CHIRON</td><td>24000</td><td>4578–8762</td><td>415</td><td>725</td></tr><tr><td>Dogsheaven*</td><td>P. Cacella</td><td>312°05′20″</td><td>0.51 m</td><td>Lhires III</td><td>5780</td><td>4550–4777</td><td>164</td><td>420</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Mean high-resolution (R ~ 80 000) CTIO-SMARTS 1.5 m/CHIRON spectrum of ζ Pup. S/N ~ 2500 in the continuum close to the vicinity of the He II λ4686 line profile.</figcaption></figure>

final homogeneous continuum normalization of the 615 spectra collected by SASER together with the 415 spectra obtained from CTIO-SMARTS 1.5 m/CHIRON and the 24 SAAO 1.9 m/GIRAFFE spectra.

characteristics of these spectra are described in section 2 of Eversberg et al. (1998).

### 2.2.4 CFHT 3.6 m/Reticon

We re-extracted the 50 archival optical spectra of ζ Pup taken with the (now decommissioned) spectrograph mounted at the coudé focus of the 3.6 m Canada-France-Hawaii Telescope (CFHT) during the nights of 1995 December 10/11 and December 12/13 (PI: Moffat) and published by Eversberg et al. (1998) in order to compare the properties of the LPVs that they observed in He II λ4686 to those that we detected during our campaign. Details on the

## 3 PHOTOMETRIC VARIABILITY

### 3.1 Amplitudes of the observed light variations

As depicted in Fig. 1, the light variations of ζ Pup appear to be coherent in the two *BRITE* filters, with qualitatively the same amplitudes. In order to quantitatively check the validity of this visual impression, we plotted the observed blue versus red amplitude variations (Fig. 3). Since the measurements from the two filters were not acquired exactly simultaneously, interpolation over the common time sampling grid of the two light curves is required in order to
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of the amplitudes of the variations of $\zeta$ Pup observed through the two *BRITE* optical filters. Top: evolution of (a) the Pearson correlation coefficient, (b) the standard deviation of the distribution of orthogonal distances and (c) the number of points involved, as a function of gap tolerance for the linear interpolation process. The vertical dashed line indicates the location of the optimal gap tolerance value (1.65 h). Middle: comparison diagrams for (left) the *BRITE* light curves presented in Fig. 1 and (right) the residual light curves after removal of the main periodic signal intrinsic to the star (Sections 3.4 and 3.5). The dashed lines indicate the reference slope = 1. Bottom: distributions of distances perpendicular to the direction of correlation. The dashed red curves are Gaussian fits to the distributions, whereas the continuous black ones are logistic distribution fits.</figcaption></figure>

construct such a diagram (which we shall hereafter call a 'BRITE-b versus BRITE-r diagram'). Therefore, a proper treatment of gaps in the time series needs to be considered, since interpolation over a large gap (such as the extreme case of the ~5-d gap during the time interval [5504–5509]) would lead to wrong estimates of the correlation coefficient, whereas if the gap tolerance is too short the number of points that remain in the diagram would be too small to allow for the extraction of meaningful slope and correlation coefficient values. The ideal case would be to consider a gap tolerance of the order of a *BRITE* orbit, which is acceptable only if the time samplings in the two light curves are such that a reasonable number of points remain for the interpolation. Thus, we scanned a grid of gap tolerance values and measured the Pearson correlation coefficient $\rho$, the standard deviation of the distribution of distances orthogonal to the direction of correlation $\sigma$, and the total number of points $N$ involved in these estimates. The first three upper panels of Fig. 3 trace the evolution of these three quantities over the considered range of gap tolerance values, and show that at a gap tolerance of 1.65 h (which is of the order of a *BRITE* orbit), $\rho$ reaches its global maximum

($\rho_{max} = 0.895$), $\sigma$ its global minimum ($\sigma_{min} = 2.23 \pm 0.05$ mmag) and $N$ is large enough to make these values meaningful ($N = 1092$). We established our diagram with this gap tolerance value (middle panels of Fig. 3). Thus, we confirm that the variations observed in the blue and the red filters are strongly correlated ($\rho = 0.895$), while the orthogonal linear regression yields a slope indicating that the amplitudes of the variations of $\zeta$ Pup in the blue filter are only $7 \pm 3$ per cent higher than the amplitudes measured through the red filter (middle left panel of Fig. 3). This relatively small difference could be due to the fact that the wavelength coverages of the *BRITE* filters both fall in the domain of validity of the Rayleigh–Jeans approximation for a hot star like $\zeta$ Pup, where different wavelengths react essentially identically to temperature changes.

Besides these considerations on the amplitudes of the light variations of $\zeta$ Pup measured in the two *BRITE* filters, it is worth noting that the distributions of orthogonal distances depicted in the bottom panels of Fig. 3 are slightly leptokurtic (excess kurtosis ~0.7), such that these orthogonal distances can be better described as logarithically distributed rather than normally distributed. More importantly, the value of the standard deviation $\sigma_{Logistic} = 2.26 \pm 0.06$ mmag is indicative of the expected scatter of instrumental origin in the *BRITE* light curves of $\zeta$ Pup. Assuming that the three satellites involved in the observations provide the same data quality, the expected scatter of instrumental origin is $\sigma_{Logistic}/\sqrt{2} = 1.60 \pm 0.04$ mmag, which is in line with the rms scatter listed in the last column of Table 2, the rms of these values being ~1.25 mmag.

## 3.2 BRITE and Coriolis/SMEI probe the stellar photosphere

In view of the typical 10–20 mmag light variations of $\zeta$ Pup as detected in the photometric observations, the next crucial point that we investigate is the stellar domain in which the observed variability originates. In principle, the variability could originate in either the photosphere or the wind of $\zeta$ Pup (or even both), since some fraction of the photospheric continuum photons interact with matter in the wind. To estimate the fraction of scattered photons in the wind, we calculate a model atmosphere for $\zeta$ Pup using the non-local thermodynamic equilibrium (non-LTE) Potsdam Wolf–Rayet (PoWR) tool, applicable for any hot star, including OB-type stars (e.g. Oskinova et al. 2011; Shenar et al. 2015; Sander et al. 2015). As described in detail by Hamann & Gräfener (2003, 2004) and Todt et al. (2015), the tool solves the radiative transfer and rate equations in expanding atmospheres under the assumption of spherical symmetry and stationarity of the flow.

To construct the model, we used the stellar and wind parameters of $\zeta$ Pup listed in Table 1: $\log(L)$, $T_{eff}$, $\log g$, $M$, $\dot{M}$, $v_{\infty}$, $\beta$ and $f_{\infty}$, all derived by Bouret et al. (2012). Fig. 4 depicts the evolution of the Rosseland mean optical depth $\tau_{Ross}$ as a function of the radial distance to the photosphere (defined at $\tau_{Ross} = 2/3$), as predicted by PoWR. We also mark the radius at which the stellar wind initiates, which is defined at the point where $v(r)$ exceeds $v_{sonic}$. It is worth noting that, in the wind domain, $\tau_{Ross}$ is virtually identical to the Thomson optical depth $\tau_{th}$ which originated in the scattering of photons off free electrons.

It is evident from Fig. 4 that the wind optical depth is $\tau_{Ross,wind} = 0.037$. This corresponds to a scattering of $1 - e^{-\tau} \approx 3.5$ per cent of the photospheric continuum light in the wind. This is therefore the maximum amount of variability that is expected to originate from scattering of photospheric photons in the wind. However, since some lines generated primarily in the wind, such as He II $\lambda$4686, are variable on the ~5 per cent–10 per cent level (Section 4.1.1, Fig. 15 and Section 4.2.1, Figs 20 and 21), the wind can
<figure><img src="image_6.png" /><figcaption>**Figure 4.** How $\tau_{\text{Ross}}$ varies as a function of the radial extent from the photosphere (defined at $\tau_{\text{Ross}} = 2/3$), as predicted by the PoWR stellar atmosphere code. Radii beyond the vertical dashed line correspond to supersonic expansion velocities, i.e. the wind domain.</figcaption></figure>

contribute at most 0.35 per cent to the variability in the continuum. This means that in the typical 10–20 mmag light variations observed by *BRITE* and *Coriolis/SMEI* in ζ Pup, only 0.035–0.07 mmag may originate from the wind. We therefore conclude that the variability we observe in the photometric measurements originates primarily in the photosphere of ζ Pup.

### 3.3 Search for periodic signals

We performed a Fourier analysis of the *BRITE* light curves of ζ Pup using the discrete Fourier transform (DFT)-based software package *PERIOD04* (Lenz & Breger 2005), suitable for the extraction of individual sinusoidal components in multiperiodic unevenly sampled time series through a pre-whitening procedure, with the possibility to perform iterative Monte Carlo simulations to assess formal uncertainties on the extracted parameters. A global view of the raw amplitude spectra of the *BRITE* light curves of ζ Pup (Fig. 5, left-hand panel) shows a gradual increase in power towards lower frequencies, suggesting the presence of a red noise component that happens to have roughly the same amplitudes in the two filters. Therefore, at each stage of the pre-whitening we computed the signal-to-noise ratio (S/N) for the detected frequency peak $v_i$ within an interval spanning [$v_i - 1; v_i + 1$] d⁻¹ which is a good compromise for a reasonable local estimate of the mean noise level while taking into account the red noise trend. Then, we only consider that the frequency peak is significant if its S/N > 4.0, a threshold that was validated both empirically (Breger et al. 1993, 1999) and from numerical simulations (fig. 4b in Kuschnig et al. 1997).

As conspicuously depicted in the amplitude spectra of the *BRITE* light curves of ζ Pup, only two significant frequency peaks are simultaneously detected in the two filters, corresponding to the 1.78-d period previously found in *Coriolis/SMEI* observations of ζ Pup, but this time its first harmonic is also prominent. This visual impression is quantified in Table 4 in which we report all the Fourier components that are either detected simultaneously in the observations through the two *BRITE* filters regardless of their statistical significance, or are significant but detected only in one filter. Only

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: DFT of the *BRITE* light curves of ζ Pup, with the insets showing zooms on the two significant frequency peaks. For ease of direct comparison, the vertical axis for the DFT of the red filter data (bottom panel) is reversed. Middle: rotation of ζ Pup as observed during the *BRITE* observing run, the value of the rotation period *P* = 1.78063(25) d assessed as the weighted mean of the values of the 1.78-d period detected during the *Coriolis/SMEI* and *BRITE* runs (Table 4). The original light curves are displayed (small grey points) as well as averages over bins of 0.04 in phase (large blue/red points) along with their $1\sigma$ uncertainties. The reference epoch HJD 2452738.7 was chosen to fall before the beginning of the *Coriolis/SMEI* 2003–2004 observing run and to allow for an easy monitoring of the evolution of the pattern of two consecutive bumps in the light curves (Section 3.4.2, Fig. 7) and the ‘S’ patterns detected in the He II λ4686 wind emission line (Section 4.1.1, Fig. 15). Right: comparison diagram for the phased light curves of ζ Pup observed through the two *BRITE* optical filters (top), together with the distribution of distances orthogonal to the direction of correlation and its Gaussian fit (bottom).</figcaption></figure>
<table><caption>**Table 4.** Fourier analysis of the *Coriolis/SMEI* and *BRITE* light curves of ζ Pup. The second column indicates the frequency resolution for each observing run. Formal 1σ uncertainties on the parameters were estimated from 100 Monte Carlo simulations. S/N values were estimated within intervals spanning 2d⁻¹ centred on each frequency. For *BRITE*, all frequencies that are detected in both the red and the blue filters are listed even if their S/N < 4.</caption><thead><tr><th>Observing run</th><th>ν<sub>res</sub><br>(d<sup>-1</sup>)</th><th>Frequency<br>(d<sup>-1</sup>)</th><th>Period<br>(d)</th><th>Semi-amplitude<br>(mmag)</th><th>S/N</th><th>Denotation</th></tr></thead><tbody><tr><td colspan="7"><em>Coriolis/SMEI</em></td></tr><tr><td>2003–2004</td><td>0.00246</td><td>0.56116 ± 0.00011</td><td>1.78202 ± 0.00035</td><td>6.64 ± 0.48</td><td>7.3</td><td>ν<sub>0</sub></td></tr><tr><td>2004–2005</td><td>0.00473</td><td>0.56148 ± 0.00019</td><td>1.78101 ± 0.00060</td><td>7.38 ± 0.50</td><td>8.3</td><td>ν<sub>0</sub></td></tr><tr><td>2005–2006</td><td>0.00569</td><td>0.56181 ± 0.00031</td><td>1.77996 ± 0.00098</td><td>6.27 ± 0.52</td><td>6.1</td><td>ν<sub>0</sub></td></tr><tr><td>2003–2006</td><td>0.00093</td><td>0.56150 ± 0.00002</td><td>1.78094 ± 0.00006</td><td>6.61 ± 0.23</td><td>12.2</td><td>ν<sub>0</sub></td></tr><tr><td colspan="7"><em>BRITE 2014–2015</em></td></tr><tr><td rowspan="10">Blue filter</td><td rowspan="10">0.00615</td><td>0.56249 ± 0.00024</td><td>1.77781 ± 0.00076</td><td>3.85 ± 0.03</td><td>8.7</td><td>ν<sub>0</sub></td></tr><tr><td>1.12087 ± 0.00102</td><td>0.89216 ± 0.00081</td><td>2.33 ± 0.16</td><td>6.7</td><td>ν'<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td>1.12694 ± 0.00120</td><td>0.88736 ± 0.00094</td><td>1.86 ± 0.17</td><td>4.1</td><td>ν''<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td>1.14385 ± 0.00100</td><td>0.87424 ± 0.00076</td><td>1.14 ± 0.03</td><td>3.4</td><td>ν<sub>1</sub></td></tr><tr><td>1.28213 ± 0.00121</td><td>0.77995 ± 0.00074</td><td>0.89 ± 0.03</td><td>2.7</td><td>ν<sub>2</sub></td></tr><tr><td>1.72701 ± 0.00236</td><td>0.57904 ± 0.00079</td><td>0.99 ± 0.03</td><td>2.9</td><td>ν<sub>3</sub></td></tr><tr><td>1.77456 ± 0.00094</td><td>0.56352 ± 0.00030</td><td>0.89 ± 0.03</td><td>2.8</td><td>ν<sub>4</sub></td></tr><tr><td>0.56250 ± 0.00017</td><td>1.77778 ± 0.00054</td><td>3.88 ± 0.22</td><td>9.6</td><td>ν<sub>0</sub></td></tr><tr><td>1.12151 ± 0.00213</td><td>0.89166 ± 0.00169</td><td>2.78 ± 0.33</td><td>7.8</td><td>ν'<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td>1.12813 ± 0.00264</td><td>0.88642 ± 0.00207</td><td>1.10 ± 0.23</td><td>3.4</td><td>ν''<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td rowspan="10">Red filter</td><td rowspan="10">0.00604</td><td>1.14242 ± 0.00087</td><td>0.87533 ± 0.00067</td><td>1.03 ± 0.25</td><td>3.0</td><td>ν<sub>1</sub></td></tr><tr><td>1.28509 ± 0.00085</td><td>0.77816 ± 0.00051</td><td>0.84 ± 0.18</td><td>2.5</td><td>ν<sub>2</sub></td></tr><tr><td>1.72947 ± 0.00170</td><td>0.57821 ± 0.00057</td><td>0.86 ± 0.20</td><td>2.6</td><td>ν<sub>3</sub></td></tr><tr><td>1.76946 ± 0.00101</td><td>0.56514 ± 0.00032</td><td>0.76 ± 0.18</td><td>2.6</td><td>ν<sub>4</sub></td></tr><tr><td>4.38282 ± 0.00768</td><td>0.22816 ± 0.00040</td><td>0.71 ± 0.20</td><td>4.3</td><td>ν<sub>5</sub></td></tr><tr><td>0.56241 ± 0.00015</td><td>1.77806 ± 0.00047</td><td>3.86 ± 0.16</td><td>10.1</td><td>ν<sub>0</sub></td></tr><tr><td>1.12111 ± 0.00047</td><td>0.89197 ± 0.00037</td><td>2.46 ± 0.36</td><td>8.0</td><td>ν'<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td>1.12714 ± 0.00078</td><td>0.88720 ± 0.00061</td><td>1.50 ± 0.34</td><td>4.3</td><td>ν''<sub>0</sub> ~ 2ν<sub>0</sub></td></tr><tr><td>1.14280 ± 0.00065</td><td>0.87504 ± 0.00050</td><td>1.05 ± 0.16</td><td>3.1</td><td>ν<sub>1</sub></td></tr><tr><td rowspan="5">Two filters<br>combined</td><td rowspan="5">0.00579</td><td>1.28323 ± 0.00678</td><td>0.77928 ± 0.00412</td><td>0.80 ± 0.15</td><td>2.6</td><td>ν<sub>2</sub></td></tr><tr><td>1.72788 ± 0.00172</td><td>0.57874 ± 0.00058</td><td>0.77 ± 0.34</td><td>2.7</td><td>ν<sub>3</sub></td></tr><tr><td>1.77440 ± 0.00090</td><td>0.56357 ± 0.00029</td><td>0.76 ± 0.17</td><td>2.7</td><td>ν<sub>4</sub></td></tr><tr><td>4.38313 ± 0.00101</td><td>0.22815 ± 0.00005</td><td>0.58 ± 0.15</td><td>4.0</td><td>ν<sub>5</sub></td></tr></tbody></table>

the 1.78-d component (denoted by $ν_0$ in Table 4) and its first harmonic have S/N > 4 and are detected in both filters. We also point out in Table 4 that a pair of frequencies, {$v'_0, v''_0$}, is closely spread around the true value of the first harmonic. It is known that a new frequency extracted at one stage of the pre-whitening procedure could be very close to a component that was found during one of the previous stages (e.g. Blomme et al. 2011). That behaviour is not surprising and is easily conceivable if for instance the physical phenomenon responsible for that signal is subject to a phase shift during the time of the observations. In all cases, one must adopt a resolution criterion in order to decide if a frequency peak is truly unique. In the present case, we use the criterion of Loumos & Deeming (1978) stating that two frequencies must be separated by at least 1.5 times the resolution frequency to be considered unique. The resolution frequency for each set of observations is also provided in Table 4. Using this criterion, we conclude that the pair {$v'_0, v''_0$} represents the same frequency, namely the first harmonic of the fundamental frequency $ν_0$.

We also note that one frequency peak, $ν_5 = 4.38313(101) d^{-1}$ [$P = 5.47560(126)$ h], is detected with a S/N = 4.3 in the observations in the red filter but undetected in the blue filter. This periodicity has never been found in previous observational campaigns on ζ Pup. At this point caution must be exercised as for the physical interpretation of this frequency. Future observations will confirm whether this signal is intrinsic to the star or an artefact of instrumental origin in our set of observations in the red filter.

Finally, we detect a set of frequencies, labelled {$ν_1; ν_2; ν_3; ν_4$} in Table 4, corresponding to periods in the range ~13.5–21 h, that are present in the two filters but with a rather low significance (S/N ~ 2.5–3.5).

Since the 1.78-d signal was discovered for the first time in *Coriolis/SMEI* observations of ζ Pup spanning 2003–2006 (Howarth & Stevens 2014), for comparison purposes and in order to explore the stability of this signal (Section 3.4.2) we also revisited the *Coriolis/SMEI* light curves and report in Table 4 the values of all frequencies detected with S/N > 4 from the Fourier analysis of the three seasonal observing runs taken separately (2003–2004, 2004–2005 and 2005–2006), and the ensemble of the observations spanning 2003–2006. We confirm that the only significant periodicity that is present in the *Coriolis/SMEI* light curves of ζ Pup is the 1.78-d period, and that the first harmonic of that period is not detected with a high significance level when taking the three seasonal observing runs separately nor when considering the entire data set as a single light curve as did Howarth & Stevens (2014). In the latter case, the values that we obtained for the fundamental frequency and its amplitude as well as their formal uncertainties are in full agreement with the values reported by Howarth & Stevens (2014) even if (1) they used a date-compensated DFT algorithm (DCDFT: Ferraz-Mello 1981) instead of the classical DFT that we adopted in our analyses, and (2) their formal 1σ uncertainties were assessed from 10 000 iterations of Monte Carlo simulations of synthetic data whereas ours were extracted from only 100 iterations.
## 3.4 The 1.78-d period

Here, we investigate the properties of the 1.78-d monoperiodic signal in ζ Pup by first determining its origin (Section 3.4.1), then tracking its evolution during the *Coriolis/SMEI* and *BRITE* observing runs (Section 3.4.2), and finally modelling the photospheric source causing it (Section 3.4.3).

### 3.4.1 Stellar pulsations versus rotational modulation

As already mentioned in Section 1, Howarth & Stevens (2014) proposed an interpretation of the 1.78-d periodic signal detected by *Coriolis/SMEI* in ζ Pup as possibly due to low-order ($0 < l \le 2$) oscillatory convection (non-adiabatic $g^−$) modes. The discovery of oscillatory convection modes with high-order azimuthal degree ($l \gtrsim 10$) in theoretical pulsation models for stars with $log L/L_⊙ = 5.0$ and $3.65 \le log T_{eff} \le 4.0$ was first reported by Shibahashi & Osaki (1981). Later the theoretical work of Saio (2011) unveiled the possibility for the excitation of low-order ($0 < l \le 2$) oscillatory convection modes in luminous stars with higher effective temperatures ($log L/L_⊙ \gtrsim 4.5$ and $log T_{eff} \gtrsim 4.0$), necessarily associated with a subsurface convection zone caused by the opacity peak due to the partial ionization of iron-group elements near $log T \sim 5.3$. In their thorough analysis of the *Coriolis/SMEI* light curves of ζ Pup, Howarth & Stevens (2014) encountered the classical dilemma ‘pulsations versus rotational modulation’ for the interpretation of the newly discovered 1.78-d signal. Such a situation is not uncommon and, depending on the amount of information available it may not be straightforward to draw a conclusion (e.g. Degroote et al. 2010, 2011; Aerts et al. 2013; Ramiaramanantsoa et al. 2014). In the case of ζ Pup, Howarth & Stevens (2014) concluded that the 1.78-d monoperiodic signal could arise from stellar oscillations rather than rotational modulation for the following reasons:

[i] the period was found to be ‘marginally consistent with the shortest possible rotation period’;

[ii] there was only a ‘minor periodogram peak’ at the first harmonic of the 1.78-d period, such that the resulting phase-folded light curve was ‘only slightly non-sinusoidal’ showing only one dominant bump (fig. 2 in Howarth & Stevens 2014);

[iii] the observed pulsation constant ($Q = P\sqrt{\bar{\rho}/\bar{\rho}_⊙} \simeq P\sqrt{M/M_⊙/(R/R_⊙)^3}$, $\bar{\rho}$ and $\bar{\rho}_⊙$ being the mean stellar and solar densities) implied by the 1.78-d period is of the order of 0.1–0.2 d, which is consistent with the predictions in the models of Saio (2011) for low-order oscillatory convection modes ($Q \sim 0.2–0.3$ d).

However, now the *BRITE* observations provide us with additional information, the most remarkable of which being the change of shape of the phased light curves (Fig. 5, middle panel; see also Fig. 7, Section 3.4.2) with respect to the single-bumped modulation observed during the *Coriolis/SMEI* observing run (fig. 2 in Howarth & Stevens 2014; see also Figs A1–A3): the *BRITE* light curves phase-folded on the 1.78-d period are highly non-sinusoidal, characterized by a pattern of two consecutive bumps separated by $\Delta\phi \sim 0.4$ which explains why the first harmonic is prominent in the Fourier spectra. The highly non-sinusoidal nature of this monoperiodic signal, coupled with the change of shape of the phased light curves during the two observing runs are strong indications that it cannot arise from pulsations but rather a signature of rotational modulation. Note that the monoperiodic and non-sinusoidal nature of the signal alone are not sufficient for drawing any conclusion on whether it comes from stellar oscillations or from rotational

modulation, as some radial pulsators are known to exhibit highly non-sinusoidal modulation of their light curves (e.g. Abt 1952, Natale, Marconi & Bono 2008, Szabó et al. 2010). Also, it is worth noting that the shape-changing property of the light curve alone could be qualified as grossly similar to the shape-changing nature of the observed light variations in RV Tau variables (e.g. Pollard et al. 1996), which are pulsators showing the period doubling effect (Moskalik & Buchler 1990). To date, all the pulsators known to exhibit this effect are radial pulsators, although in theory the pulsation modes involved do not necessarily need to be radial but need to satisfy the half-integer resonance criterion, i.e. two modes A and B satisfying $(2n + 1)v_A = 2v_B$ with *n* integer (Moskalik & Buchler 1990). In that case, due to resonant phase locking, mode B is not observed as an independent frequency: in the amplitude spectrum, mode A manifests at its frequency $v_A$ along with its harmonics, whereas mode B manifests through subharmonic frequencies of $v_A$. Therefore, if the 1.78-d variability in ζ Pup were to be explained in terms of the period doubling effect, the main pulsation mode A has to be the one oscillating at the ~21.4-h period, while mode B would manifest as the 1.78-d period. But then, in that situation, the fact that there are no signs of any subharmonics of the ~21.4-h period (except the 1.78-d period) in any of the observing runs remains unexplained, and more importantly the ~21.4-h period itself is quasi-absent during the ~4 yr of *Coriolis/SMEI* observing run (fig. 2 in Howarth & Stevens 2014; Figs A1–A3). Based on all these considerations, we conclude that the highly non-sinusoidal nature of this monoperiodic shape-changing 1.78-d signal in ζ Pup is incompatible with stellar oscillations but can be naturally explained as arising from rotational modulation due to the presence of bright spots appearing and disappearing at different locations on the stellar surface.

The pattern of two consecutive bumps observed in the phase domain during the *BRITE* run is noteworthy. Patterns of two consecutive dips are common in the non-sinusoidally modulated light curves of cool stars having two dominant dark surface spots (e.g. McQuillan, Aigrain & Mazeh 2013, Davenport, Hebb & Hawley 2015), whereas patterns of two consecutive bumps are common in the phased light curves of stars exhibiting bright surface chemical inhomogeneities (Bruntt et al. 2009; Bernhard et al. 2015a; Bernhard, Hümmerich & Paunzen 2015b; Weiss et al. 2016). Over the past decade, an increasing number of light curves of massive OB stars have been seen to show low-frequency non-sinusoidal signal best explained as arising from rotational modulation rather than pulsations as the corresponding phased light curves show a pattern of two consecutive bumps and the periodograms show prominent peaks at the harmonics of the fundamental frequency (De Cat & Aerts 2002; Degroote et al. 2010, 2011; Balona 2016). In the case of ζ Pup, one noticeable property in the frequency domain is the prominence of the fundamental frequency (*Coriolis/SMEI* and *BRITE* observing runs) and the first harmonic (*BRITE* observing run), while higher harmonics are not detected. This behaviour is reminiscent of the findings of Balona (2016) who scanned a sample of OB stars observed during the K2 mission and found that 54 per cent of them show rotational modulation rather than pulsations: 39 per cent of these rotational variables show a single dominant low-frequency peak in their periodogram (no harmonics) and 61 per cent have both the fundamental frequency and its first harmonic prominent in their periodogram while higher harmonics are absent. This particular behaviour reconciles with the theoretical investigations of Clarke (2003) who found that the Fourier series decomposition of stellar light variations due to surface spots is characterized by significant contributions from the fundamental and the first harmonic, with the
<table><caption>**Table 5.** Estimates of the minimum possible rotation period $P_{\text{min}}$, the equatorial velocity $v_e$ and the inclination angle $i$ of $\zeta$ Pup, taking $P = 1.78063(25)$ d for the value of the rotation period and considering the two cases where $R_e = \sqrt{1.5}R$ and $R$. The values of $P_{\text{min}}$ were calculated for $g_p = g$.</caption><thead><tr><th>Parameter</th><th></th><th>$R_e = \sqrt{1.5}R$</th><th>$R_e = R$</th></tr></thead><tbody><tr><td>$P_{\text{min}}$</td><td>(d)</td><td>$2.10^{+0.48}_{-0.43}$</td><td>$1.90^{+0.43}_{-0.39}$</td></tr><tr><td>$v_e$</td><td>(km s<sup>-1</sup>)</td><td>$66.1 \pm 132$</td><td>$539 \pm 108$</td></tr><tr><td>$i$</td><td>(°)</td><td>$19.4^{+7.3}_{-4.7}$</td><td>$24.0^{+9.4}_{-5.9}$</td></tr></tbody></table>

strengths of the higher harmonics rapidly vanishing as the order increases.

Under these considerations, now that our analyses indicate that the 1.78 signal in $\zeta$ Pup is related to rotational modulation rather than pulsations, we revisit the three points [i]–[iii] mentioned above. We note that since the calculation of the pulsation constant mentioned in point [iii] involves the stellar mass and radius, which have large uncertainties in the case of $\zeta$ Pup due to the large uncertainty in its distance, the last point [iii] can only be confirmed or ruled out once a refinement on the value of the distance is available, e.g. with the Gaia mission (Lindegren et al. 2012; Martín-Fleitas et al. 2014). Also, as previously discussed in detail, point [ii] is an integral part of the key property that the shape of the light curve changes due to the appearance of spots at different locations on the stellar surface. Finally, point [i] also depends on the stellar radius which is affected by the uncertainty in the distance. As pointed out by Howarth & Stevens (2014), in order to maintain a positive equatorial effective gravity in a Roche model, the stellar rotation period must not be shorter than

$$P_{\text{min}} = 3\pi \sqrt{\frac{R_e}{g_p}}, \qquad (1)$$

where $R_e$ is the equatorial radius and $g_p$ the gravity at the pole. Note that to date the available values of $R$ and $g$, which are listed in Table 1, were all derived from spherical models (Martins, Schaerer & Hillier 2005; Bouret et al. 2012). The values of $P_{\text{min}}$ for $g_p = g$ and considering the two cases where $R_e = \sqrt{1.5}R$ and $R_e = R$ are summarized in Table 5, indicating that the 1.78-d period is close to the minimum possible value of the stellar rotation period, but not entirely excluded by virtue of the estimated uncertainties. Also, $\zeta$ Pup being a fast rotator, if effects of gravity darkening are non-negligible, the value for $P_{\text{min}}$ could be lower as the polar gravity would be higher than the value of the effective surface gravity $\log g(\text{cm s}^{-2}) = 3.64 \pm 0.1$ listed in Table 1. This is clearly illustrated on Fig. 6: for instance in the case $R_e = \sqrt{1.5}R$, if $\log g_p$ is 10 per cent higher than $\log g = 3.64$, then the minimum rotation period drops to $1.38^{+0.34}_{-0.29}$ d. From all these considerations, we conclude that the three points [i]–[iii] are compatible with the fact that the 1.78-d signal in $\zeta$ Pup comes from rotational modulation.

Given the values of the 1.78-d period yielded by the Fourier analyses of the different data sets of $\zeta$ Pup (Table 4), we adopt a unique value of the rotation period throughout our investigation by assessing the weighted average of the values of the 1.78-d period detected in the three Coriolis/SMEI seasonal observing runs and the combined BRITE light curve. We use the inverse square of the $1\sigma$ uncertainties derived from the Monte Carlo simulations as weights on the periods, yielding a final weighted average of $P = 1.78063(25)$ d. Such a relatively short rotation period for $\zeta$ Pup is quite surprising at first sight, but is more consistent with the suspected rotational evolution of the star involving past interactions within a massive bi-

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Minimum rotation period of $\zeta$ Pup (equation 1) as a function of the dimensionless fraction $x = (\log g_p/3.64) - 1$ which represents the percentage of deviation of the polar gravity from the value of the effective surface gravity $\log g(\text{cm s}^{-2}) = 3.64 \pm 0.1$ derived by Bouret et al. (2012). Two cases are considered: $R_e = \sqrt{1.5}R$ (solid blue) and $R_e = R$ (dashed red). In each case, the shaded area delimits the upper and lower bounds imposed by the uncertainties on the stellar radius and a 2.75 per cent uncertainty on $\log g$. The dashed-dotted horizontal line is the period $P = 1.78063(25)$ d (its $1\sigma$ uncertainty is of the order of the thickness of the line itself).</figcaption></figure>

nary (Vanbeveren, van Rensbergen & De Loore 1998; see also Section 5 and Fig. 22) or even a multiple system (Vanbeveren 2012). In all cases, the relatively short rotation period has strong implication on both the value of the stellar inclination angle $i$ with respect to the line of sight and the equatorial velocity $v_e$. Given the available estimates of the stellar radius $R = 18.99 \pm 3.80 R_⊙$ and the projected rotational velocity $v_e \sin i = 219 \pm 18 \text{ km s}^{-1}$ (Table 1), the stellar inclination angle can be directly assessed as:

$$i = \arcsin \left( \frac{P v_e \sin i}{2\pi R_e} \right), \qquad (2)$$

while the equatorial velocity is:

$$v_e = \frac{2\pi R_e}{P}. \qquad (3)$$

The values of these two parameters are also reported in Table 5 for the two cases $R_e = \sqrt{1.5}R$ and $R$. As clearly depicted by the numbers in Table 1, we can conclude that the short rotation period implies that the star is seen at a very low inclination angle, but bearing in mind that all these values strongly depend on the stellar radius which remains highly uncertain as a result of the large uncertainty in the distance to the star.

### 3.4.2 Evolution of the 1.78-d signal

As noted in the previous section, the key property leading to the conclusion that the 1.78-d signal in $\zeta$ Pup is related to rotational modulation rather than stellar oscillations is the change of shape of the phased light curve from the epoch of the Coriolis/SMEI observations to the epoch of the BRITE run. This behaviour led us to also investigate the evolution of the amplitude of the 1.78-d signal in $\zeta$ Pup through the BRITE observing run and through each Coriolis/SMEI seasonal run. In order to reduce the effects of high-frequency stochastic variability intrinsic to the star (that we shall characterize in Section 3.5), we applied an infinite impulse response low-pass 32nd-order Butterworth filter to the light curves,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Investigation of the evolution of the 1.78-d signal in ζ Pup during the *BRITE* observing run (2014–2015). Time increases upwards. Left: time–frequency diagram of the combined filtered *BRITE* light curves constructed from a 12-d sliding window Fourier transform. The upper panel shows the power spectrum of the entire unfiltered light curve (red line) compared to that of the entire filtered light curve (black line), along with the magnitude response of the low-pass filter (dashed line) and the location of its cut-off frequency (vertical dotted line). Middle: continuous dynamic phased light curves of ζ Pup. The upper panel shows the phased light curve for the entire unfiltered data (black filled circles) compared to that from the entire filtered light curve (green filled circles). Right: average shapes of the phased light curve of ζ Pup (black diamonds) over the first 20 rotational cycles (Part I, see the text), the next 10 rotational cycles at the beginning of the transition phase (Part II), the last 15 rotational cycles during the transition phase (Part III), the next 25 rotational cycles (Part IV) and the last 27 rotational cycles (Part V). These shapes are compared to the overall shape of the phased light curve (green filled circles).</figcaption></figure>

cutting off at $v_c = 1.4 \text{ d}^{-1}$. It is known that Butterworth filters have the advantage of having a monotonic magnitude response that is optimally flat (no ripples) in the passband, which is the prime reason why we chose it over other types of filters, at the cost of a less steep roll-off that is mitigated by increasing the order of the filter. We inspected the reliability of our low-pass Butterworth filter in studying the 1.78-d signal in ζ Pup by comparing the power spectra of the original and the filtered light curves (Fig. 7, upper left panel). There is globally no significant difference between the two power spectra in the filter bandpass, the loss of power being at the ~2 per cent level at the fundamental frequency peak $v_0$ and ~8 per cent at the first harmonic $2v_0$. Besides, in the time domain, the rms scatter of the residuals between the phased light curve before and after applying the low-pass filter is 0.26 mmag, which is of the order of half the rms value of the $1\sigma$ uncertainties of the mean points in the original phased light curve (Fig. 7, upper middle panel). From these results, we conclude that the low-pass filtering process is reliable to study the evolution of the 1.78-d signal, while reducing the contribution of the high-frequency stochastic variations intrinsic to the star which shall be discussed in detail in Section 3.5. Thus, we performed all of the analyses related to the 1.78-d signal on the filtered light curves.

In the frequency domain, we conducted a Fourier analysis of the *BRITE* light curve in sliding windows over restricted time intervals spanning $\Delta T = 12 \text{ d}$ (with a step of $\Delta T \times 1$ per cent), thus covering at least ~6 stellar rotational cycles at each step. The lower left panel of Fig. 7 depicts the resulting time-frequency diagram, clearly

showing that the fundamental frequency is always present from the beginning of the observing run until almost the end (around day 5610.0), while the first harmonic only emerges from day ~5510.0 and remains visible until the end of the observing run. Between days ~5540.0 and ~5575.0, the fundamental frequency is still present but its power is partly absorbed by the first harmonic which is very prominent during that time interval.

Besides the time-frequency analysis, we also tracked the evolution of the 1.78-d signal by splitting the light curve into consecutive rotational cycles and stacking them on top of each other to form the grey-scale diagram in the middle panel of Fig. 7 as a function of rotational phase and cycle number. Such a diagram is not common practice for analysing periodic/cyclic signals in photometric observations of stars as it requires a large number of points per cycle, but could be useful in detecting and characterizing subtle variations (e.g. Gies et al. 2012, Davenport et al. 2015). In our case, the combined filtered *BRITE* light curve of ζ Pup provides us with an increased number of points per rotational cycle and increased S/N so that the diagram turns out to be a very useful complement to the time-frequency analysis as it helps us track events in the light curve that explain the behaviour that we found in the time-frequency diagram. More specifically, before day ~5490 the phased light curves show one large bump centred around phase 0.6, which is the reason why we only see the fundamental frequency and not the first harmonic during that epoch of the observations. This behaviour could be due to the presence of one relatively large bright spot or a group of smaller bright spots at the photosphere. Then the time between
days ~5490 and ~5538 is characterized by a transition that eventually leads to the formation of two bumps separated by $\Delta\phi \sim 0.44$ in the phased light curve at the end of the transition. These two bumps are clearly defined between days ~5540.0 and ~5575.0, causing the prominence of the first harmonic of the 1.78-d signal during that time interval in the time-frequency diagram. The presence of these two bumps could be the manifestation of two dominant bright spots at the stellar surface, separated by 360° $\times \Delta\phi \simeq 158^\circ$ in longitude. The separation of these two bumps in phase seems to become larger towards the end of the observing run, a sign that the surface spots experience a slight migration in longitude, which is then an appropriate explanation for the detection of the pair of frequencies $\{v'_0; v''_0\}$ closely spread around the true value of the first harmonic during the pre-whitening procedure as pointed out in Section 3.3. The right-hand panel of Fig. 7 shows the mean phased light curves over several consecutive rotational cycles, thus summarizing the above descriptions:

[i] the first 20 rotational cycles observed before the transition phase, showing one large bump centred around phase 0.6 (Cycles 2396–2415, which we will hereafter refer to as Part I),

[ii] the next first 10 rotational cycles during the transition phase (Cycles 2416–2425; Part II),

[iii] the last 15 rotational cycles during the transition phase (Cycles 2426–2440; Part III),

[iv] the 25 rotational cycles right after the transition phase, during which the two bumps separated by $\Delta\phi \sim 0.44$ are clearly visible (Cycles 2441–2465; Part IV),

[v] the last 27 rotational cycles of the observing run (Cycles 2466–2492; Part V).

We note that the behaviour of the star during the transition (Parts II and III) is not clearly defined at this point, the two most pertinent scenarios being either that the large spot/spot group in Part I integrates and the resulting most dominant parts migrate to form the two dominant spots that we see in Part IV at locations on the surface separated by ~158° in longitude, or the large spot/spot group in Part I gradually disappears while two new dominant spots take birth.

Howarth & Stevens (2014) already investigated the stability of the 1.78-d signal during the *Coriolis/SMEI* observing run by performing a DCDFT on the seasonal subsets of the light curve and on 50-d subsets with 50 per cent overlaps. They noticed that the fundamental frequency of the signal does not change, but its amplitude varies by a factor 2 during the observing run, while it experiences a modest phase excursion (fig. 3 in Howarth & Stevens 2014). In order to investigate the stability of the 1.78-d signal on the 12 yr time-scale covered by the *Coriolis/SMEI* and the *BRITE* observing runs, we established the same dynamic plots as Fig. 7 for the *Coriolis/SMEI* seasonal observing runs (Figs A1–A3). The time-frequency plots show that the fundamental frequency is always detected except at some points where it is so weak that its power is comparable with the noise level (e.g. between HJD−2451545 ~ 1450−1500). Furthermore, the dynamic phased light curves show that the minima wander within the interval 0.4–0.6, confirming the behaviour reported by Howarth & Stevens (2014). More interestingly, the phased light variations observed during the epoch of the *Coriolis/SMEI* observing run present only one relatively large bump located around phase zero, which could indicate that the surface spots seen during the *Coriolis/SMEI* observing run appear at different longitudes compared to those observed during the epoch of the *BRITE* observations. However, at this point, a strict interpretation of the position of the maxima in the phased light curve cannot be ascertained, since the exact value of the rotation period remains unknown, and

the value of its error bar $\sigma = 0.00025$ d could lead to an apparent phase shift of ~0.03 in 1 yr, or ~0.35 over the ~12 yr covered by the *Coriolis/SMEI* and the *BRITE* observing runs. Nevertheless, even considering such a large phase shift, the clear shape-changing nature of the light variations over these 12 yr, characterized by the presence of one large bump during the *Coriolis/SMEI* observing run and the presence of two dominating bumps during the epoch of the *BRITE* observations, is already a strong indication that the spots seen in 2003–2006 are not the same as the ones seen in 2014–2015.

### 3.4.3 Mapping the spotted stellar surface

Approaches to the modelling of light curves showing effects of rotational modulation associated with the presence of surface spots can be classified into two distinct categories depending on whether the problem is solved in a direct manner or in an inverse manner. The direct method consists in calculating the emerging stellar light variations produced by a given number of spots (often taken as a fixed parameter), provided some a priori assumptions on their shape (often assumed to be circular), their decay law, their size and their location as well as the time of their first appearance on the stellar surface. Then these last three free parameters are refined for each spot until the calculated light curve reasonably fits the observed light curve. Obviously, in that method a convergence towards a reasonable fit to the observed light curve strongly depends on the assumed number of spots. In practice, one could gradually increase the assumed number of spots until the smallest number that produces a good fit is reached, and then argue that in some sense the simplest acceptable spot distribution has been found. But then the other five parameters for each of the spots are still being assumed. This direct method, also known as the analytical approach, has been used in the modelling of light curves of cool low-mass stars exhibiting active surface regions (e.g. Kipping 2012; Béky et al. 2014; Davenport et al. 2015) and a few cases of hotter stars (Lüftinger et al. 2010; Degroote et al. 2011; Ramiaramanantsoa et al. 2014), offering the advantage of being fast in terms of computation time but suffering from the disadvantage that the number of spots, their shapes and their decay law have to be assumed a priori along with good estimates of starting points for the other three free parameters for each spot. Alternatively, the inverse approach (Harmon & Crews 2000; Roettenbacher et al. 2011, 2013, 2016; Krtička et al. 2007; Kolbin, Shimansky & Sakhibullin 2013; Kolbin & Shimansky 2014; Kolbin, Sakhibullin & Gabdeev 2015) consists in assessing the distribution of specific intensity at the stellar surface that reasonably fits the observed light curve, with no a priori assumptions concerning the six parameters mentioned previously, the sole assumption being that the stellar surface contains spot-like features. This conceptually constitutes an inverse ill-posed problem and has been proven to require more computation time than the direct approach (thus its designation as the ‘numerical approach’), keeping in mind its great advantage of being free from any prior assumption on the five spot parameters enumerated previously.

In order to model the bright spot-induced light variations of $\zeta$ Pup as observed by *BRITE* and *Coriolis/SMEI*, we adopted the numerical approach by using a constrained non-linear light-curve inversion algorithm (LI:<fnref n="2" /> Harmon & Crews 2000). This algorithm

<footnote n="2">² Initially denominated ‘matrix light-curve inversion’ in Harmon & Crews (2000) as it was inspired from a matrix-based algorithm for the inversion of light curves of planetary objects to assess their surface albedo distribution (Wild 1989, 1991). However, the actual implementation of the algorithm no</footnote>
for mapping the surface of spotted stars was initially formulated and tested on simulated data by Harmon & Crews (2000), and applied to the analyses of the evolution of dark spots on the K2IV primary component of II Pegasi (Roettenbacher et al. 2011) and cool K-type Kepler targets (KIC 5110407, Roettenbacher et al. 2013; KOI-1003, Roettenbacher et al. 2016). Here, we apply the algorithm for the first time to map the locations and evolution of bright spots at the surface of a hot massive star.

**Light-curve inversion – the algorithm.** Let us represent the observed light curve as the time series $I = \{I(t_k)\}_{1 \le k \le K}$, thus containing K measurements. The algorithm as it stands assumes rigid rotation, and thus currently does not explicitly allow for surface differential rotation. Then the basis of LI is to partition the stellar surface into N latitude bands of equal angular size ($\Delta\vartheta = \pi/N$), each nth latitude band containing $P_n$ spherical rectangles having exactly equal area. In addition, the partitioning is made such that all the patches on the stellar surface are of nearly equal area, which is achieved by having $P_n$ proportional to the cosine of the latitude of the centre of the band, subject to the constraint that $P_n$ must be an integer. Thus, at a given time $t_k$, if say $J_{np}$ is the contribution of the specific intensity along the outward normal for patch (n; p) that subtends a solid angle $\Omega_{np}(t_k)$ as seen from the observer, and $\mathcal{L}_{np}(t_k)$ the value of the limb darkening at the location of that patch, then to the limit of a large total number of patches, the intensity from the star at that time takes the discretized form

$$I(t_k) = \sum_{n=1}^{N} \sum_{p=1}^{P_n} \Omega_{np}(t_k) \mathcal{L}_{np}(t_k) J_{np}, \quad (4)$$

in which the outer summation needs only to be performed up to $n_s \le N$, which is the index of the southernmost latitude band that is visible by the observer as a consequence of the inclination of the star (thus $n_s = N$ if the star is seen equator on). Then the idea is to find the set of relative patch intensities $\hat{J} = \{\hat{J}_{np}\}_{1 \le n \le n_s, 1 \le p \le P_n}$ that produces a reconstructed light curve $\hat{I} = \{\hat{I}(t_k)\}_{1 \le k \le K}$ such that the rms deviation (expressed in magnitudes) between $\hat{I}$ and the observed light curve $\hat{I}$ is equal to the estimated rms noise in $\hat{I}$ (also expressed in magnitudes). The ill-posed nature of the problem clearly appears here, the key component that is involved being the noise in the observations that shows as high-frequency ripples that would be indistinguishable from the effect of a multitude of small bright and dark spots randomly distributed over the stellar surface. Simply trying to find the distribution of patch intensities that best fits the observed light curve, i.e. through a standard minimization of a goodness-of-fit criterion, would inevitably yield a granulated surface dominated by noise artefacts. As in any ill-posed problem, a regularization process must be involved in order to mitigate noise sensitivity and ensure convergence towards a unique solution (Twomey 1977; Craig & Brown 1986). In our case, regularization is achieved by solving for the distribution of patch intensities $\hat{J}$ through a Lagrange constrained minimization of an objective function defined as a Lagrangian of the form

$$E(\hat{J}, \mathbf{I}, \lambda, B) = G(\hat{J}, \mathbf{I}) + \lambda S(\hat{J}, B), \quad (5)$$

in which $G(\hat{J}, \mathbf{I})$ is the term that measures the goodness of fit between the reconstructed and the observed light curves, while the second term $\lambda S(\hat{J}, B)$ controls the smoothness of the reconstructed stellar surface. The function $S(\hat{J}, B)$ is therefore called smoothing

longer involves matrices; thus, we adopt its common current designation: LI.

function, such that the lower value it takes the smoother the surface is. Thus, in this context, the Lagrange multiplier $\lambda$ acts as a smoothing parameter, a property that can be easily understood when considering the two asymptotic cases where $\lambda \to 0$ and $+\infty$: the first case amounts to minimizing $E(\hat{J}, \mathbf{I}, \lambda, B) \simeq G(\hat{J}, \mathbf{I})$ which would yield solutions dominated by noise artefacts as discussed above, while the second case would yield very smooth solutions that poorly fit to the observed light curve as the term $\lambda S(\hat{J}, B)$ would dominate over $G(\hat{J}, \mathbf{I})$. The optimal solution of the inverse problem corresponds to an intermediate value of $\lambda$ that is the best trade-off between these two limit cases, yielding a reconstructed light curve $\hat{I}$ that reasonably fits the observed light curve $\mathbf{I}$, while mitigating the tendency of overfitting.

The subtle role of the smoothing function $S(\hat{J}, B)$ can be understood from its closed-form expression which is a generalized Tikhonov regularizer:

$$S(\hat{J}, B) = \sum_{n=1}^{n_s} \sum_{p=1}^{P_n} w_n c_{np} [\hat{J}_{np} - \langle \hat{J} \rangle]^2, \quad (6)$$

where $\langle \hat{J} \rangle$ is the average patch intensity, $w_n$ are latitude-dependent weighting factors and the coefficients $c_{np}$ incorporate the constraint that the surface exhibits confined discernible spots on a uniform background. This is achieved through the definition:

$$c_{np} = \begin{cases} 1 & \text{if } \hat{J}_{np} < \langle \hat{J} \rangle \\ B & \text{if } \hat{J}_{np} \ge \langle \hat{J} \rangle \end{cases}, \quad (7)$$

with B a strictly positive real constant. Keeping in mind that $S(\hat{J}, B)$ is involved in the objective function $E(\hat{J}, \mathbf{I}, \lambda, B)$ to be minimized, two distinct cases arise from this definition of $c_{np}$: $B > 1$ and $0 < B < 1$. The first case would mean that a patch that tends to be brighter than $\langle \hat{J} \rangle$ by a given amount rather than being dimmer than $\langle \hat{J} \rangle$ by the same amount will increase $S(\hat{J}, B)$ by a factor of $B > 1$, thus inflicting a penalty in the minimization of $E(\hat{J}, \mathbf{I}, \lambda, B)$. Conversely, if $0 < B < 1$, the penalty becomes a gain in the minimization of $E(\hat{J}, \mathbf{I}, \lambda, B)$ as the contribution of the patch to $S(\hat{J}, B)$ decreases by a factor of $1/B > 1$. In other words, the information on whether we want to model dark or bright spots on a uniform background is incorporated in the coefficients $c_{np}$: the modelling of dark spots requires $B > 1$ as was adopted in the work of Roettenbacher et al. (2011, 2013, 2016), while in this investigation, we constrain $0 < B < 1$ to model the bright spots on $\zeta$ Pup. Note that for these reasons, $S(\hat{J}, B)$ is also called the penalty function and B the bias parameter, since it biases the solution in favour of either dark or bright spots.

Also, the presence of the latitude-dependent weighting factors $w_n$ is noteworthy and comes from the necessity to take into account the fact that the maximum projected area of a patch increases the closer to the sub-Earth latitude it is located. Hence, the amount of light modulation induced by any given spot is equivalent to that of a smaller spot located closer to the sub-Earth latitude. This degeneracy would lead the algorithm to a systematic convergence towards solutions with spots lying near the sub-Earth point, as a result of the fact that small spots incur small values of the penalty function $S(\hat{J}, B)$. To reduce that systematic behaviour, the nth latitude band is weighted by a coefficient $w_n$ that is proportional to the difference between the maximum and minimum values of the product of the projected area and the limb darkening for patches in the band, namely:

$$w_n = \frac{1}{W} \times \left| \max_{1 \le p \le P_n} \{A_{np}\} - \min_{1 \le p \le P_n} \{A_{np}\} \right|, \quad (8)$$
in which $W = \max_{1 \le n \le n_s} \{w_n\}$ is a normalization factor, while $A_{np}$ is the product of the projected area and the limb-darkening value for the patch $(n; p)$:

$$A_{np} = \frac{1}{2} \left[ \sin i (\cos \varphi_1 - \cos \varphi_2) \left( \frac{\pi}{N} + \frac{1}{2} (\sin 2\varphi_1 - \sin 2\varphi_2) \right) + \frac{2\pi}{P_n} \cos i (\cos^2 \varphi_1 - \cos^2 \varphi_2) \right] \times \mathcal{L}_{np}, \quad (9)$$

in which $i$ is the stellar spin axis inclination angle with respect to the line of sight, $\varphi_1$ and $\varphi_2$ are the colatitudes of the northern and southern edges of the latitude band, whereas $\varphi_1$ and $\varphi_2$ are the longitudes of the western and eastern edges of the patch. As a matter of fact, $\max_{1 \le p \le p_n} \{A_{np}\}$ is reached at the sub-Earth meridian, while $\min_{1 \le p \le p_n} \{A_{np}\}$ is reached at the anti-Earth meridian. Although it was originally found that the introduction of the latitude weighting coefficients $w_n$ generally resulted in a modest improvement on the latitude resolution and an increased computation time (Harmon & Crews 2000), here we decided to keep them in order to optimize the output of LI.

Finally, the goodness-of-fit term $G(\hat{\mathbf{J}}, \mathbf{I})$ is defined to be the rms deviation in magnitudes of the reconstructed light curve from the observed light curve. Evaluated in units of the estimated noise variance in magnitudes $\sigma^2$ in the observed light curve, the closed-form expression for $G(\hat{\mathbf{J}}, \mathbf{I})$ is:

$$G(\hat{\mathbf{J}}, \mathbf{I}) = \frac{(2.5 \log e)^2}{K \sigma^2} \sum_{k=1}^{K} \left( \frac{I(t_k) - \hat{I}(t_k)}{I(t_k)} \right)^2. \quad (10)$$

As mentioned earlier, the goal of LI is to find the optimal distribution of patch intensities $\hat{\mathbf{J}}$ that yield a reconstructed light curve $\hat{\mathbf{I}}$ such that the rms variance between $\hat{\mathbf{I}}$ and the observed light curve $\mathbf{I}$ reaches $\sigma^2$. In other words, an inversion consists in finding the optimal values of the Lagrange multiplier $\lambda$ and the bias parameter $B$, for which we have:

$$G(\hat{\mathbf{J}}(\lambda; B), \mathbf{I}) = 1. \quad (11)$$

That is achieved through a two-step process involving a Van Wijngaarden-Dekker-Brent root-finding algorithm (Press et al. 2007) to first solve equation (11) for $\lambda$ for a set of values of the bias parameter $B$ constrained within the interval $[0; 1]$, yielding $\lambda_{\text{opt}}(B)$ which is in turn used to find $B_{\text{opt}}$ as the root of the function

$$g(B) = \frac{\min\{\hat{J}_{np}(\lambda_{\text{opt}}(B); B)\}}{\langle \hat{J}_{np}(\lambda_{\text{opt}}(B); B) \rangle} - \frac{J^{(s)}}{J^{(p)}}, \quad (12)$$

where $J^{(s)}$ and $J^{(p)}$ are, respectively, the spot and photosphere specific intensities, two parameters that we determine by evaluating the Planck function at the central wavelength of the filter passband for the spot and photosphere temperatures. In principle, a better way to obtain a proxy for these two parameters is to integrate the Planck function over the filter passband for the spot and photosphere temperatures, the most ideal way being a numerical integration of the stellar spectral energy distribution over the filter passband. However, given the typical uncertainty in spot temperatures, the method that we adopted is enough for achieving a nearly optimal precision.

**Spot temperatures.** At this point, particular consideration has to be given to spot temperatures: since our BRITE observations were performed in two different filters, any amplitude difference between observations through the two filters could be exploited in order to extract potential information on the typical temperatures of the

spots. This can be easily conceived by assuming to first order that all the spots that are present on the surface at a given time $t_k$ are all at the same temperature $T^{(s)}$ (a reasonable approximation if the spots roughly probe the same stellar atmosphere layer), and assuming that their combined area at that snapshot in time is equivalent to a single spot occupying a fractional area $f$ on the stellar surface subtending a limb angle $\theta$ as seen from the centre of the star and emitting $I^{(s)}$ per unit area, while the unperturbed stellar photosphere emits $I^*$ per unit area in that filter. This situation amounts to the configuration considered by Clarke (2003) (see his fig. 1), so that the observed photometric signal takes exactly the same closed-form expression as his equation (1):

$$F_j(\theta) = \mathcal{K} \left[ I_j^* + f \left( I_j^{(s)} - I_j^* \right) \cos \theta \right], \quad (13)$$

in which we have the additional subscript $j$ to indicate observations through filter $j$ ($j$ = 'b' or 'r' in our case). The factor $\cos \theta$ takes into account the foreshortening of the projected area of the spot according to the angle $\theta$, whereas the factor $\mathcal{K}$ is a constant that depends on the stellar radius and the distance to Earth. Furthermore, equation (13) can be expressed in terms of magnitudes and approximated to first order as:

$$m_j - m_j^* \approx -2.5(\log e) \left( \frac{I_j^{(s)} - I_j^*}{I_j^*} \right) f \cos \theta. \quad (14)$$

with $m_j^* = -2.5 \log (\mathcal{K} I_j^*)$ the stellar magnitude in the absence of spots. Therefore, the ratio of the observed variations in the blue filter to the observed variations in the red filter is:

$$\alpha = \frac{\Delta m_b}{\Delta m_r} = \frac{m_b - m_b^*}{m_r - m_r^*} \approx \frac{I_b^* - I_b^{(s)}}{I_r^* - I_r^{(s)}} \times \frac{I_r^*}{I_b^*}. \quad (15)$$

But by definition, this ratio $\alpha$ turns out to be exactly the slope in our BRITE-b versus BRITE-r diagrams! In terms of colour indices, equation (15) yields:

$$(b-r)_d = (b-r)_* - 2.5 \log \alpha, \quad (16)$$

where $(b-r)_d$ is the colour index of the difference between the spot and the photosphere, while $(b-r)_*$ is the colour index for the unperturbed star. In other words, equation (16) remarkably tells us that a simple measurement of the slope in the BRITE-b versus BRITE-r diagram can give us an estimate of spot temperatures by means of calibrations between $(b-r)$ and effective temperatures. Also, equation (16) can be applied to any dual-band photometric time series observations of any star when appropriate calibrations of colour indices versus effective temperatures are available. Unfortunately, in our situation, we measure a slope $\alpha = 0.99 \pm 0.09$ in our BRITE-b versus BRITE-r diagram for the 1.78-d signal in $\zeta$ Pup (right-hand panel of Fig. 5), which does not allow us to extract any information on spot temperatures. This insensitivity could be due to the fact that the two BRITE red and blue filters fall in the domain of validity of the Rayleigh-Jeans approximation for a hot star like $\zeta$ Pup. Thus, given the fact that there is no amplitude difference in the observations in the two filters, we performed our inversions on the combined BRITE light curve of $\zeta$ Pup, allowing us to have decreased errors in the mean points and reduced gaps in rotational phase coverage, even if it has been demonstrated by Harmon & Crews (2000) that the inversion of light curves in multiple filters can improve the spot latitude resolution.

**Mapping the surface of $\zeta$ Pup during the BRITE observing run.** Under all the considerations in the subheaders of Section 3.4.3, we performed inversions of the BRITE light curves of $\zeta$ Pup by partitioning the stellar surface into $N = 60$ latitude bands with the
first four equatorial bands containing 90 patches, and we adopt as input parameters:

[i] the estimated rms noise σ for the observed light curves split into different parts (Section 3.4.2). As mentioned earlier, the goal of LI is to find the ‘smoothest’ solution until the rms residual between the reconstructed and observed light curves reaches the level of σ. It has been demonstrated that in practice (Harmon & Crews 2000; Roettenbacher et al. 2011, 2013, 2016), the reconstructed surface will show obvious noise artefacts when the assumed noise level becomes too low. Our criterion for choosing the best solution is thus defined by the effective noise level at which this behaviour starts to occur. Therefore, we only assessed the rms noise in our observed light curves in order to choose the appropriate range of assumed values of σ that need to be considered. We measured σ in the range 1.2–0.7 mmag in the phased light curves of ζ Pup observed by BRITE, hence we scanned the following range of σ in our inversions: 1.5, 1.2, 0.9, 0.8, 0.7, 0.6, 0.5 and 0.4 mmag;

[ii] the inclination angle: we performed inversions for $i = 24°$, which is the mean value defined by the limits that we found in Section 3.4.1 (Table 5), and we also performed inversions for the lower limit $i_{min} = 15°$ and the upper limit $i_{max} = 33°$;

[iii] the photosphere temperature (which we take to be $T_{phot} = T_{eff} = 40.0$ kK: Table 1), and a range of discrete values for the spot temperature: $T^{(s)} = 41.5, 42.5, 45, 47.5$ and $50$ kK;

[iv] limb-darkening coefficients: we adopted a quadratic limb-darkening law and extracted the appropriate coefficients for the stellar parameters of ζ Pup in the Bessel V filter (the closest to the middle of BRITE combined filters) from the grid of limb-darkening coefficients calculated by Reeve & Howarth (2016) using non-LTE, line-blanketed TLUSTY model atmospheres.

The first important point that we noticed was that the inversions for the assumed range of spot temperatures only slightly differ in the resulting spot-to-photosphere contrast: higher spot temperatures resulted in slightly brighter spots (which is expected for the reasonable values of spot temperature that we considered), but negligibly affected the locations and the shapes of the spots found by the algorithm. Figs 8–10 illustrate the outcome of our inversions of the BRITE light curve of ζ Pup for $T^{(s)} = 42.5$ kK and for $i = 24°$, $i = i_{min} = 15°$ and $i = i_{max} = 33°$, respectively. We clearly note that the range within which the stellar inclination angle varies does not influence the locations of the detected spots on the stellar surface. However, we notice that, as the inclination angle gets lower, the shapes of the spots tend to be more stretched in longitude. In other words, the spots are rounder for $i = 33°$ than for $i = 15°$. The equatorial spots appear to be the most affected by this behaviour, which could then be explained as follows. Given a certain value of the stellar inclination angle, consider an equatorial spot having a certain shape that provides the right amount of light modulation at a given rotational phase: if the stellar inclination is decreased, the lower part of the spot would become less visible (its overall projected area on to the plane of the sky decreases), such that the spot has no other choice than to flatten in order to compensate for that behaviour and provide the same amount of light modulation as in the case of the higher inclination angle. Also, in general it is expected that there should be a tendency for the spots to shift to higher latitudes as the assumed inclination angle decreases, since the sub-Earth latitude moves northward in that case. However, for small equatorial spots, which is the case for most of the spots detected by the algorithm here, that tendency will not be as pronounced because a small equatorial spot will be visible for half the rotation period regardless of the value of the inclination angle. This explains why

the range within which the stellar inclination angle varies here does not affect the locations of the spots detected by the algorithm. In all cases, our inversions of the BRITE light curve of ζ Pup converge towards the following behaviour:

[i] Part I: one spot $S_1$ is detected lying close to the equator, along with a slightly larger spot $S_2$ located at mid-latitude (~55°:5). Their separation is ~66° in longitude.

[ii] Part II: the two spots $S_1$ and $S_2$ seem to be subject to a shift in longitude, while a third spot ($S_3$) appears close to the equator, separated from the previous equatorial spot $S_1$ by ~158° in longitude.

[iii] Part III: spot $S_2$ has almost faded away, while the equatorial spots $S_1$ and $S_3$ remain there. A fourth spot $S_4$ seems to be located mid-way between $S_1$ and the previous location of $S_2$. We identify three possible scenarios to explain the presence of $S_4$: either it is exactly $S_2$ that is migrating towards $S_1$, or part of the disintegration of $S_2$ migrating towards $S_1$, or a completely new spot.

[iv] Part IV: only two dominant equatorial spots, $S_1$ and $S_3$, are visible separated by ~158° in longitude.

[v] Part V: the two dominant equatorial spots at $S_1$ and $S_3$ remain visible. Both spots seem to have slightly shifted in longitude (by ~6° for $S_1$ and ~16° for $S_3$), and their separation becomes ~148° in longitude.

Thus, our inversion of the BRITE light curve of ζ Pup has enabled us to map the surface of the star, find the locations of the surface bright spots that best explain the observed variations, and track the evolution of these surface inhomogeneities. We also performed inversions of the light curves of ζ Pup obtained during the three Coriolis/SMEI seasonal runs. The resulting surface maps for $T^{(s)} = 42.5$ kK and $i = 24°$ are illustrated on Figs B1–B3. We detect on average ~2–3 spots at the equator and at mid-latitudes during the epoch of the Coriolis/SMEI observations. The spots clearly appear at different longitudes compared to the BRITE observing run, but whenever there are two equatorial spots, their separation seems to be ~148°–180° in longitude.

## 3.5 Stochastic variability

As mentioned in Section 3.3, the sole significant period that emerged from our Fourier analysis of the BRITE and Coriolis/SMEI light curves of ζ Pup was the 1.78-d signal due to rotational modulation (with the first harmonic of the fundamental period also prominent in the amplitude spectrum of the BRITE light curves). In order to check for any signs of other types of variability in the light curves, we therefore removed this 1.78-d periodic signal. Since the shape of this signal changes in time, the standard removal by pre-whitening is the least reliable method as it assumes a constant amplitude and phase. Instead of calculating the residual light curve by pre-whitening the 1.78-d signal, we performed its removal in the phase domain by subtracting templates of the phased light curves from the different parts of the observing runs: in the case of the BRITE observations, we adopted the five different parts of the observing run (Parts I–V; left-hand panel of Fig. 7), while in case of the Coriolis/SMEI observations, we adopted the different parts defined on the left-hand panels of Figs A1–A3. Then we recomputed the amplitude spectra of the residual light curves to confirm that there is no sign of the 1.78-d signal left and no new significant periodicity introduced during the removal process, thus confirming our previous findings in Section 3.3. Now a close inspection of the residual BRITE light curves in the two filters reveals variations reaching ~20 mmag in
<figure><img src="image_10.png" /><figcaption>**Figure 8.** LI: mapping the photosphere of ζ Pup as observed by BRITE in 2014–2015, for $T^{(s)} = 42.5$ kK and a stellar inclination angle $i = 24^\circ$. Time increases upwards. The left-hand panel illustrates the observed light curve (filled circles) during Parts I–V of the BRITE observing run as defined on Fig. 7, Section 3.4.2, along with the reconstructed light curve (green line), with the residuals plotted below the light curves. Then follows a view of the star at rotational phase 0.375 (middle panel) and the pseudo-Mercator projection of the stellar surface (right-hand panel). The vertical open brackets on the left of the pseudo-Mercator projections indicate the range of latitudes visible by the observer. The sub-Earth point is at longitude 0° at rotational phase zero.</figcaption></figure>

peak-to-peak amplitudes, stochastically generated, but coherent for several hours (Fig. 11). The first hint indicating that this stochastic signal is intrinsic to the star is the fact that it behaves the same way in the observations in the two *BRITE* filters as clearly seen in Fig. 11. This is further confirmed by the behaviour of the *BRITE-b* versus *BRITE-r* diagram for the residual light curves in which we measure a strong correlation ($\rho = 0.856$) and a slope $\alpha = 1.10 \pm 0.04$ (right-hand panel of Fig. 3). Also as already shown in Fig. 11, our ultimate check that confirms that this stochastic signal is not of instrumental origin nor due to an artefact introduced during the decorrelation process is a comparison with the *BRITE* light curve of *ι Car* (HD 80404, A8Ib, V=2.25), which is of similar magnitude as ζ Pup, was also observed by *BRITE* as part of the Vela/Pup observing run, has data that were decorrelated with the same decorrelation

routines, and shows no obvious signs of periodic variability from our Fourier analysis. From all these considerations, we conclude that this stochastic signal is intrinsic to the star. We emphasize here the fact that the stochastic nature of this signal describes the way each feature in the signal appears to be generated, while a given feature seems to be well organized (or coherent) for several hours after its generation. Also from these considerations, we can take the liberty to call this signal ‘noise’, as long as we keep in mind that it is a noise intrinsic to the star and not of instrumental origin. It is thus worth noting that this stochastic behaviour is reminiscent of the findings of Balona (1992) who came to the conclusion that ζ Pup is an irregular microvariable (with amplitudes ~20 mmag) from ground-based Strömgren *b* filter photometric monitoring of the star.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as in Fig. 8 but for *i* = *i*<sub>min</sub> = 15°. The locations of the detected spots remain the same as for *i* = 24° (Fig. 8), while their shapes appear to be more stretched in longitude.</figcaption></figure>

The origin of this stochastic signal remains unknown, with poten-
tial candidates being unresolved randomly excited p-mode oscilla-
tions similar to those discovered in the O-type primary component
of the massive binary system HD 46149 (Degroote et al. 2010), or
IGWs generated from either a subsurface convection zone (Cantiello
et al. 2009) or from the convective core (Aerts & Rogers 2015), the
latter case having been proposed to be to date the best interpretation
of the red noise components in the CoRoT observations of the three
O-type dwarfs HD 46223, HD 46150 and HD 46966 in the young
open cluster NGC 2244 (Blomme et al. 2011). Also, perhaps the
set of frequencies {ν₁; ν₂ν₃; ν₄} (P ~ 13.5–21 h) detected by our
Fourier analysis in both the observations in the red and blue filters
but with low significance (S/N ~ 2.5–3.5; Section 3.3, Table 4) is
part of this stochastic signal, which would explain their low signif-
icance as they would be present/excited only occasionally and last
for a limited time. In all cases, regardless of its origin, in view of its

typical amplitudes the quantification of this stochastic variability in
terms of its amplitudes and time-scales is as important as the char-
acterization of the 1.78-d periodic rotationally induced modulation
itself.

3.5.1 Stochastic variability – amplitudes

One striking property of this stochastic variability in ζ Pup is that
it reaches ~20 mmag in peak-to-peak amplitude during the epoch
of the BRITE observations (Fig. 11 and right-hand panel of Fig. 3),
while we should recall that the periodic 1.78-d signal does not
exceed ~15 mmag peak to peak (Fig. 5). The overall standard
deviation in the amplitudes of the residual BRITE light curves of ζ
Pup is σ_brite ≃ 6 mmag. But as assessed in Section 3.1, the scatter
of instrumental origin is σ_brite, i = 1.60 ± 0.04 mmag. Since that
instrumental scatter is much smaller than the overall scatter in the
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as in Fig. 8, but for *i* = *i*max = 33°. The locations of the detected spots remain the same as for *i* = 24° (Fig. 8), while their shape appear to be rounder.</figcaption></figure>

residual light curve, we can infer that the standard deviation of the amplitudes associated with the stochastic signal intrinsic to the star during the epoch of the *BRITE* observations is of the order of 6 mmag.

Regarding the *Coriolis/SMEI* observations of ζ Pup, the overall standard deviation in the amplitudes of the residual *Coriolis/SMEI* light curve of the star is higher: σ*SMEI* ≃ 16.8 mmag. However, the contribution of instrumental effects to this scatter remains ambiguous. Indeed, in view of the *Coriolis/SMEI* observations of the late O-type subgiant ζ Oph (O9.2IVnn; *V* = 2.6) having an rms dispersion of ~20 mmag dominated by instrumental effects, while contemporaneous *MOST* observations of the star has a dispersion of ~6 mmag, most of which is intrinsic to the star (Howarth et al. 2014), it is impossible to constrain the nature of the stochastic component of the intrinsic variability of ζ Pup during the epoch of the *Coriolis/SMEI* observations.

### 3.5.2 Stochastic variability – time-scales

As mentioned earlier, the features constituting the intrinsic stochastic signal in ζ Pup seem to be coherent for several hours. In order to characterize the time-scales of these transient features, given the random nature of their generation we performed a continuous wavelet analysis of the residual combined *BRITE* light curves. We tested two cases: (1) application of the 1.10 factor on the amplitudes of the observations in the red filter prior to the generation of the combined light curves, and (2) no amplitude correction prior to the combination of the light curves. We noticed no significant difference in the resulting continuous wavelet transforms (CWT), thus we chose to work on the combined light curves without any amplitude correction. Also it is important to note that in general, the continuous formalism can only be applied on well sampled and equally spaced data. Our overall combined *BRITE* light curve of ζ Pup is well
<figure><img src="image_13.png" /><figcaption>**Figure 11.** A 15-d subset of the residual *BRITE* (2014–2015) light curve of ζ Pup after removal of the 1.78-d periodicity related to rotational modulation, along with the simultaneously observed *BRITE* (2014–2015) red light curve of ι Car (HD 80404, A8Ib, V=2.25) used as a comparison star.</figcaption></figure>

sampled enough, but obviously not regularly sampled. However, by choosing appropriate subsets of the combined light curve that are densely sampled and resampling these subsets with one single sampling interval, it is reasonable to adopt the continuous formalism. Therefore, our approach is to subdivide the combined light curve into subsets that contain no gaps greater than 5 h (39 subsets), then refining the selection by keeping only subsets spanning at least 2.5 d, so that we ended up with 11 subsets with time bases between ~3–52 d and mean sampling intervals in the range 53.2–75.4 min. Lastly, through a shape-preserving piecewise cubic interpolation (Sprague 1990), we resampled all the 11 subsets with one single sampling interval taken to be the minimum of the mean sampling intervals of the subsets ($\delta t_s$ = 53.2326 min). Once the subsets are defined, we come to the choice of an appropriate wavelet family that will help us optimally and robustly characterize the time-scales of the transient features composing the stochastic signal intrinsic to ζ Pup. Given the facts that the observed variations in the light curve seem to equally have bumps and dips between ±10 mmag (Figs 11 and 3), and also since we want to have the most reliable conversion from wavelet scales to pseudo-periods, the wavelet family that is the most appropriate to our case is the family of real-valued Morlet wavelets, for which the mother wavelet is a sinusoid modulated by a Gaussian envelope. Furthermore, this definition itself clearly indicates that this family of wavelets is the most suitable one for the characterization of sinusoidal signals with finite lifetimes, which could be the case of the observed stochastic variations here if it is related to unresolved stochastically excited pulsations with finite lifetimes.

Fig. 12 shows an example of the outcome of our CWT of a 30-d subset of the combined residual *BRITE* light curve of ζ Pup, indicating a gradual increase in power towards longer time-scales, followed by a drop in power for time-scales longer than ~17 h. The bump at a pseudo-period of ~16.5 h might correspond to the 16.67 ± 0.81 and 16.90 ± 0.48 h periods reported by Berghoefer et al. (1996) from an 11 d simultaneous X-ray and H$\alpha$ monitoring of ζ Pup. These ‘periods’ have never been found again since then, and Nazé et al. (2013) invoked the possible transient nature of the signal associated with them. In our present investigations, the

scalogram on Fig. 12 clearly shows that the ~16.5 h pseudo-period is only sporadically present during the *BRITE* observing run, most prominently around time ~5519.5 and ~5535.0, and lasts about ~2.5–4 d when it is excited. Also given the width of the bump at the ~16.5 h pseudo-period, it is not relevant to assign an exact value for it, as for instance a Gaussian fit to the bump yields 16.42 ± 9.5 h. In view of the scalogram on Fig. 12, we notice a pseudo-period drift from long time-scales (~22 h) to shorter time-scales (~12 h), mostly visible between times ~5517.5 and 5521.0 during which there is a pseudo-period drift from ~16 to ~11 h, and between times ~5533.0 and 5536.0 during which pseudo-periods drift from ~22 h to ~16 h, a behaviour that remains unexplained but could be the cause for the large spread of power around the bump at the ~16.5 h pseudo-period. Finally, we also note a minor pseudo-period bump at ~8 h, which also seems to be only intermittently excited during the observing run, but has a shorter lifetime (≤ 1 d) than the ~16.5 h pseudo-period. From all these considerations, it could also be that these two bumps are just part of the spectrum of pseudo-periods constituting the stochastic variability of the star. Therefore, from the CWT of all the 11 subsets that we analysed, we calculated the weighted mean energy distribution of the time-scales (with the lengths of the subsets as weights), best fit with a power law of index 3.17 ± 0.06 for time-scales less than 8 h, while time-scales greater than 8 h follow a less steep power law of index 0.37 ± 0.06 (Fig. 13).

# 4 SPECTROSCOPIC VARIABILITY

## 4.1 Large-scale wind structures

### 4.1.1 Searching for periodicity in the He II λ4686 line

As a first step in our investigations of the variability in the He II λ4686 wind emission line of ζ Pup, we searched for periodicity in the global time series of difference profiles of the line (with respect to the global average profile) by performing a CLEAN Fourier analysis (Roberts, Lehar & Dreher 1987) for each Doppler velocity in the range ±1275 km s⁻¹ (Fig. 14). The resulting power spectrum unveils a dominant peak at $v_{0,\text{He II}} = 0.5615(1) \text{ d}^{-1}$ [$P_{0,\text{He II}} = 1.7809(3)$
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Wavelet analysis of a 30-d subset of the combined residual *BRITE* light curve of ζ Pup. The main panel depicts the scalogram (percentage of energy for each wavelet coefficient) obtained from a Morlet wavelet family-based CWT of the analysed 30-d portion of light curve illustrated on the top panel. The right-hand panel shows the sum of the percentages of energy per wavelet coefficient over time, thus representing the overall energy distribution of pseudo-periods detected during this subset of the observations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Weighted average energy distribution of pseudo-periods resulting from the CWT of the 11 subsets in the combined *BRITE* light curve of ζ Pup (taking as weights the time bases of the subsets), showing a change of slope at a pseudo-period of ~8 h. The continuous lines are linear fits for pseudo-periods less than 8 h (green) and pseudo-periods beyond 8 h (red).</figcaption></figure>

d], for which we assessed the uncertainty by extracting the HWM of a Gaussian fit to the peak. This frequency peak, detected across ±500 km s⁻¹ in Doppler velocity and mostly prominent on the blue side of the line profile, corresponds to the rotation frequency of the star that was also detected in the *BRITE* and *Coriolis/SMEI* observations (Section 3.4; Table 4 and Figs 5 and 7). Moreover, the first harmonic $2v_0$.He II is visible in the power spectrum but with a lower strength. Also some of the aliases of these two peaks are still visible in the power spectrum, such as the peak at $|2v_0$.He II - 2| \sim 0.877$ d⁻¹ as well as the 1 d⁻¹ alias peak, all due to the unavoidable sampling gaps from the ground-based observations. Apart from the

rotation period and its first harmonic we do not detect any other significant periodicity in the He II λ4686 line during our campaign. The absence of any sign of the 5.1–5.2-d period through the entire width of the line profile is noteworthy. Furthermore, we do not detect any signs of that period in H α which was only covered by our CTIO-SMARTS 1.5 m/CHIRON observations. That period was previously detected in H α (Moffat & Michaud 1981) and in UV observations (Howarth et al. 1995). However, in view of our findings that the true stellar rotation period is the 1.78-d period and not the 5.1-5.2-d period (Section 3.4), the previous interpretations of the latter needs to be revised.

In view of the detection of the stellar rotation period in the variations of the He II λ4686 line profile, we inspected the behaviour of the difference spectra phase-folded on the rotation period, adopting the value *P* = 1.78063(25) d that we derived from the *Coriolis/SMEI* and *BRITE* photometric observations (Section 3.4.1). As we know that the bright spots inducing this signal in the light curve evolve during the observing campaign (Sections 3.4.2 and 3.4.3), we subdivided the spectroscopic observations according to the five different parts of the *BRITE* observing run, Parts I–V, and looked at the phase-folded dynamic difference spectra for each of these parts. Due to insufficient phase sampling, the phased dynamic difference spectra during the first part of the transition (Part II) and towards the end of the campaign (Part V) did not show any obvious pattern that could be used to characterize the nature of the LPVs related to the stellar rotation in the He II λ4686 line. The phased dynamic spectra for Part I (before the transition), Part III (end of the transition) and Part IV are illustrated on the left-hand panel of Fig. 15 in which we note the following behaviour:

[i] Part I: an ‘S’ pattern in excess emission is mildly visible, reaching maximum redshift around phase 0.35 and maximum blueshift around phase 0.85.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Global time series analysis of the He II λ4686 wind emission-line profile of ζ Pup. Top: the 1054 profiles of the line as recorded throughout the entire campaign (grey) along with their average (black). Middle: difference profiles with respect to the average profile. Main: CLEAN periodogram per Doppler velocity within ±1275 km s⁻¹ for the time series of difference spectra, obtained after 10 CLEAN iterations. The side panel depicts the normalized integrated periodogram power within Doppler velocities ±500 km s⁻¹ where most of the LPVs are detected.</figcaption></figure>

[ii] Part III: the behaviour of the LPVs are fuzzy, but the clearest one is the bump on the red side of the line profile (between 0 and 400 km s⁻¹) before phase ~0.5, which appears to be part of the 'S' pattern seen in Part I.

[iii] Part IV: a double 'S' bumped pattern is strikingly visible, wandering between ~±400 km s⁻¹. Note that this corresponds to the part during which the two bumps in the BRITE light curves are clearly visible, caused by the two dominant equatorial surface spots separated by ~158° in longitude.

The double 'S' pattern that we clearly see during Part IV is a typical sign of the presence of two CIR arms induced by two equatorial perturbations (Dessart & Chesneau 2002), which, according to our findings in Section 3.4.3, are the equatorial bright spots separated by ~158° in longitude observed by BRITE during this part of the observing campaign (Fig. 8). This situation confirms the theoretical

predictions of Dessart & Chesneau (2002), as not only did they use the stellar parameters of ζ Pup in all their CIR models, but more importantly our findings for the behaviour in Part IV are remarkably similar to the behaviour that they found in their Model C (fig. 3 in Dessart & Chesneau 2002) for which the assumed photospheric perturbations were two equatorial bright spots diametrically opposed on the stellar surface and the assumed stellar inclination angle was 20° (the actual inclination angle of ζ Pup is i ~ 24° according to our findings: Section 3.4.1), with the assumed line emission region (LER) extending from 2R to 4R.

### 4.1.2 CIR model

Following the detection of the patterns of manifestations of CIR arms in the phased dynamic difference profiles of the He II λ4686 line, most prominently the double 'S' pattern that we see during Part IV of the observing run, we further investigated the consistency of these observed LPVs with models of CIR arms that would be generated by the bright photospheric spots observed by BRITE and mapped by LI. To this end, we adopted an analytical approach to determine the behaviour of LPVs induced by a given number of CIR arms in a hot-star wind emission-line profile. In our approach, the input parameters consist of:

[i] the inclination angle i of the star, its rotation period P and its radius R;

[ii] the radius $r_I$ at which the stellar wind initiates, as well as the terminal wind speed $v_∞$;

[iii] two radial functions, one for localizing the LER, and another one for accounting for ionization gradients through the wind;

[iv] the number of CIR arms, each of them specified by its density contrast $\eta$ with respect to the unperturbed wind, its opening angle $\psi$ and the location of its photospheric driver on the stellar surface.

Here, we provide a brief description of the basis of this CIR modelling approach and the results that we obtained in the case of ζ Pup, but a more detailed description of the code will be available in a dedicated forthcoming paper (St-Louis, Tremblay & Ignace 2017).

**The emergent line profile for a spherically symmetric wind.** In the absence of CIR arms, our calculations of hot-star wind emission-line profiles is performed under the assumption of spherical symmetry, stationarity and smoothness of the outflow, the latter meaning a stellar wind free of small-scale inhomogeneities. Also in terms of the velocity field of the wind material, the polar and azimuthal components are assumed to be negligible compared to the radial expansion. Thus, with the geometry depicted in Fig. 16, a portion of wind material located at point P(r, $\vartheta$, $\phi$) in terms of its spherical coordinates in the stellar reference frame moves outward with a velocity $\vec{v}(\vec{r}) = v_t(r)\vec{e}_r$, assumed to follow the usual $\beta$-law:

$$v_t(r) = v_\infty \left(1 - b \frac{r_1}{r}\right)^\beta, \quad b = 1 - \left(\frac{v_I}{v_\infty}\right)^{\frac{1}{\beta}} \quad (17)$$

with $v_I$ the radial speed at which the stellar wind initiates. As mentioned in Section 3.2, $v_I = v_{sonic}$ and is reached at $r_I \simeq 1.027R$ (Fig. 4). Thus, the wind domain being in the supersonic regime, calculations of line profiles originating from the wind can be performed under the Sobolev approximation: at a given frequency displacement $\nu$ with respect to the line-centre frequency $\nu_0$, the emission is dominated by the contribution from locations in the stellar wind where the local component of the velocity field is in resonance with the frequency $\nu$ (see e.g. equations 38 and 39 in
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Arms of CIRs in the wind of ζ Pup as traced in its He II λ4686 emission-line profile in parallel with the BRITE 2014–2015 observing run. Time increases upwards. Left: observed dynamic difference spectra of ζ Pup phase-folded on the rotation period *P* = 1.78063(25) d, during Part I (bottom), Part III (second from the bottom) and Part IV (third from the bottom), along with the corresponding light variations observed by BRITE illustrated on the side panels. Small horizontal ticks on the right-hand side of each dynamic diagram indicate the phase sampling. The average profiles of the He II λ4686 line during each of these three parts of the observing campaign are depicted in the top panel (dashed: Part I; dotted: Part III and solid: Part IV). Right: models of CIRs corresponding to each of the three parts of the observing campaign, with an assumed stellar inclination angle *i* = 24° and taking the locations of the bright spots detected by the LI algorithm (Section 3.4.3; Fig. 8) as input locations of the photospheric footprints of the CIRs. The top panel depicts the average modelled profiles corresponding to the three parts of the observing run (dashed: Part I; dotted: Part III and solid: Part IV), not accounting for the resonance zone effect as explained in Section 4.1.2.</figcaption></figure>

Petrenz & Puls 1996). At a given frequency *ν*, the monochromatic line flux as received by the observer at distance *D* from the star is:

$$F_{\nu} = \frac{1}{D^2} \int_0^{2\pi} \int_0^{\pi} \int_0^{\infty} j_{\nu}(r) r^2 dr d\theta d\phi, \qquad (18)$$

involving the emissivity *j*<sub>ν</sub>, which, under the Sobolev approximation, is expressed as (Rybicki & Hummer 1983):

$$j_{\nu}(r) = \mathcal{F}(r)\mathcal{H}\left(\frac{v_{\nu}(r)}{c}v_{0} - |v - v_{0}|\right), \qquad (19)$$

where the presence of the Heaviside function *H* accounts for the resonance zones mentioned previously. This also reflects the fact

that at a given radial distance *r*, the spherical shell of radius *r* expanding at velocity *v*<sub>ν</sub>(*r*) contributes to the emission as a flat-top profile spanning *v*<sub>0</sub>[1 ± *v*<sub>ν</sub>(*r*)/*c*], weighted by the quantity:

$$\mathcal{F}(r) = \frac{1}{2v_0} \frac{c}{v_\nu(r)} \kappa(r)S(r)g(r)\mathcal{P}(r), \qquad (20)$$

the factor *c*/2*v*<sub>0</sub>*v*<sub>ν</sub>(*r*) arising from wind broadening, *κ*(*r*) being the integrated line opacity, *S*(*r*) the source function, *g*(*r*) a function that incorporates all dependencies on temperature variations and changes in ionization throughout the wind, while *P*(*r*) denotes the Sobolev escape probability.

For recombination lines such as the He II λ4686 line in ζ Pup, the product *κ*(*r*)*S*(*r*) scales as the square of the density *ρ*(*r*) of the wind.
<figure><img src="image_18.png" /><figcaption>Figure 16. Geometry adopted for the simulation of hot-star wind emission-line profiles from a spherically symmetric smooth stationary stellar wind and calculation of LPVs induced by the presence of a spiral-shaped CIR compression.</figcaption></figure>

By means of the equation of mass continuity, $\dot{M} = 4\pi r^2 v_r(r) \rho(r)$, the product $\kappa(r)S(r)$ therefore ultimately scales as $1/r^2 v_r^2(r)$. As for the function $g(r)$, we assume a dependence of the form $[1 - v_r(r)]^{\alpha_0}$, where the exponent $\alpha_0$ would then counter the $1/v_r^2(r)$ dependence to create an asymmetric bell curve in velocity space that localizes the LER. From the investigations of Hillier et al. (2012), we know that one has to go at least up to ~2.5–3.5R to account for the actual shape of the He II $\lambda$4686 line profile in $\zeta$ Pup, and the emission falls off slowly beyond that bulk of the wind (Hillier, private communication). Here, we adjust the exponent $\alpha_0$ to roughly localize the LER, although we note that in reality the departure coefficients from conditions of LTE may change rapidly and may not be well represented by the simple form of the $g(r)$ function that we adopted. However, it is a reasonable approximation to adopt in the absence of any detailed constraints on the LER.

In addition, the Sobolev escape probability $P(r)$ is given by:

$$P(r) = \frac{1 - e^{-\tau(\vec{r}, \vec{n})}}{\tau(\vec{r}, \vec{n})}, \qquad (21)$$

$\tau(\vec{r}, \vec{n})$ being the optical depth evaluated at the point of consideration along the direction specified by the unit vector $\vec{n}$, expressed in the Sobolev approximation as:

$$\tau(\vec{r}, \vec{n}) = G(r) \times \frac{\rho^2(r)}{|\nabla_{\vec{n}}(\vec{v} \cdot \vec{n})|}, \qquad (22)$$

where $|\nabla_{\vec{n}}(\vec{v} \cdot \vec{n})|$ denotes the velocity gradient along direction $\vec{n}$, while the function $G(r)$ accounts for ionization gradients through the wind. Following Lamers, Cerruti-Sola & Perinotto (1987), we adopt the generalized form $G(r) = v_r^{\alpha_1}(r)[1 - v_r(r)]^{\alpha_2}$, in which the exponents $\alpha_1$ and $\alpha_2$ are in practice adjusted to match the observed shape of the line profile. Thus, taking also into account the equation of mass continuity, the optical depth takes the form:

$$\tau(\vec{r}, \vec{n}) = M \times \frac{v_r^{\alpha_1 - 2}(r) [1 - v_r(r)]^{\alpha_2}}{r^4 |\nabla_{\vec{n}}(\vec{v} \cdot \vec{n})|}, \qquad (23)$$

with $M = (\dot{M}/4\pi)^2$. Then, in order to evaluate the Sobolev optical depth, the velocity gradient is measured along the line of sight ($\vec{n} = \vec{e}_z$ here) and takes the form (e.g. Rybicki & Hummer 1983; Lamers et al. 1987; Mazzali 1990; Petrenz & Puls 1996; Ignace & Hendry 2000; Ignace & Gayley 2003):

$$\nabla_{\vec{n}}(\vec{v} \cdot \vec{n}) = \frac{\partial v_z}{\partial z} = \frac{\partial v_r(r)}{\partial r} \cos^2 \theta + \frac{v_r(r)}{r} \sin^2 \theta, \qquad (24)$$

with $\theta$ denoting the observer polar angle (i.e. the angle between the line of sight, which is the z-axis here, and the position vector $\vec{r}$; Fig. 16).

Finally, we note that the calculation of the line flux from equations (18) and (19) has to be performed excluding part of the stellar wind material located behind the star which remains unseen by the observer. Occultation happens if the angle $\theta$ is such that:

$$\cos \theta < -\sqrt{1 - \left(\frac{R}{r}\right)^2}, \qquad (25)$$

in which $\cos \theta$ can be expressed in terms of the azimuth angle $\varphi$, the colatitude $\vartheta$ and the stellar inclination angle $i$ through the cosine rule:

$$\cos \theta = \cos \vartheta \cos i + \sin \vartheta \sin i \cos \varphi. \qquad (26)$$

**Presence of CIR arms.** In our simulations, the stellar surface is assumed to rotate rigidly at a constant angular frequency $\Omega = 2\pi/P$, and the photospheric source of a CIR arm is assumed to be a bright surface perturbation having a circular shape located at $S(R, \vartheta_s, \varphi_s + \Omega t)$. The half-opening angle $\psi$ of the CIR is defined to be that of the cone formed by the circular spot and the centre of the star (Fig. 16). As already pointed out by Dessart & Chesneau (2002), the feature in the CIR that is at the origin of the variability observed in the wind emission lines is the spiral-shaped region of density compression associated with the CIR (region III in fig. 5 of Cranmer & Owocki 1996), which in turn can be perceived as a region where the emissivity is increased by a factor $(\eta + 1)^2$, with $\eta = (n_{CIR} - n_{sph})/n_{sph}$ being the density contrast between the region of CIR compression of density $n_{CIR}$ and the spherically symmetric smooth wind of density $n_{sph}$. Thus, within the region of CIR compression, the integrand in the calculation of the line flux (equation 18) becomes $r^2(\eta + 1)^2 j_v$. Also, in terms of its shape, the region of CIR compression is taken to follow a spiral with increasing radial distance from the centre of the star, while at a given radial distance $r$ the intersection of the spiral-shaped CIR compression and a spherical shell of radius $r$ is assumed to form a circle. The axis of the region of CIR compression is then defined to be the locus composed by the centres of these circular cross-sections with increasing radial distance $r$. From all these considerations and from the geometry described in Fig. 16, the condition for falling within the region of CIR compression is straightforward:

$$\arccos \chi < \psi, \qquad (27)$$

where $\chi$ is the angle between the position vector $\vec{r}$ and the outward normal at the centre M of the circular cross-section of the CIR compression at radial distance $r$ (unit vector $\vec{\mu}$ at point M in Fig. 16). Again, by means of the cosine rule, this angle can be expressed in terms of the colatitude $\vartheta$ and the azimuth angle $\varphi$ of the position vector, as well as the colatitude $\vartheta_{CIR}$ and the azimuth angle $\varphi_{CIR}$ of the centre of the cross-section of the CIR arm at distance $r$ from the centre of the star:

$$\cos \chi = \cos \vartheta \cos \vartheta_{CIR} + \sin \vartheta \sin \vartheta_{CIR} \cos (\varphi_{CIR} - \varphi), \qquad (28)$$
in which $\vartheta_{CIR} = \vartheta_s$ at all times and radial distances, whereas the azimuth angle $\varphi_{CIR}$ carries all the information on the time dependence (owing to the stellar rotation) and the shape of the CIR arm:

$$ \varphi_{CIR} = \varphi_s + \Omega t - \frac{\sin \vartheta_s}{x_p} \left( \frac{r}{R} - 1 + b \ln \left[ \frac{r}{R} \right] + \frac{1-b^2}{b} \ln \left[ (1-b) \frac{v_\infty}{v_r(r)} \right] \right), \quad (29) $$

a closed-form expression that follows Ignace, Hubrig & Schöller (2009) and Ignace, St-Louis & Proulx-Giraldeau (2015) in which it is assumed that $v_\psi(r) = 0$ for the purpose of solving the spiral shape, due to the complexity of non-radial force considerations and because we do not perform hydrodynamical calculations in our simulations and do not have any real guidance on $v_\psi(r)$. Also, it is worth noting that in equation (29), the dimensionless quantity $x_p = v_\infty/(\Omega R)$ measures the ratio of the distance travelled by the outer boundary of the stellar wind in one rotation period with respect to the stellar equatorial diameter, thus indicative of the asymptotic pitch angle of the spiral-shaped CIR compression.

**Application to the case of ζ Pup.** First and foremost, we emphasize again that we generate all our CIR models to check the consistency of the observed CIR patterns (left-hand panel of Fig. 15) with the bright photospheric spots observed by BRITE and mapped by LI. No hydrodynamical calculations and no formal fit of the modelled LPVs to the observed ones are performed here. Moreover, the models have infinite S/N, as no noise of instrumental origin or intrinsic to the star (e.g. due to clumps) was added to the models.

Regarding the observed shape of the He II λ4686 emission-line profile in ζ Pup, the presence of a notch close to the line centre is noteworthy. Since the line is mainly formed in the stellar wind, the radial expansion of the wind itself can contribute to a blueshifted P Cygni-type absorption in the line. In that respect, the shape of the line profile can be considered as of type III in Beals' (1953) classification of P Cygni-type line profiles. Additionally, ζ Pup being a fast rotator, it has been shown that a proper consideration of the azimuthal velocity field can account for the actual shape of some of its photospheric emission-line profiles as well as its wind emission lines such as Hα and He II λ4686 (Hillier et al. 2012), the main mechanism involved being the so-called 'resonance zone effect' as described in detail by Petrenz & Puls (1996), who found that the introduction of the azimuthal velocity field induces a decrease of the optical depth for Doppler velocities near the line centre as the corresponding resonance zones are twisted away from the star, whereas the optical depth is increased for higher Doppler velocities as the corresponding resonance zones are twisted towards the star (fig. 6 in Petrenz & Puls 1996). As a consequence, the line emission flux is redistributed towards higher Doppler velocities, resulting in the actual shape of the He II λ4686 line having two maxima near ±v sin i and a (slightly blueshifted) notch close to the line centre. Therefore, that notch is a feature associated with the unperturbed wind, and not related to the presence of variable small- or large-scale wind structures. Thus, in our simulations we do not take into account the resonance zone effect and consider a velocity field dominated by radial expansion as already stated in Section 4.1.2.

The right-hand panel of Fig. 15 depicts the modelled patterns of LPVs corresponding to CIR arms that would be induced by the bright spots detected by the LI algorithm for Parts I, III and IV

of the observing campaign (Section 3.4.3, Fig. 8). These models correspond to a stellar inclination angle of 24° and the stellar wind parameters of ζ Pup listed in Table 1, the values of the density contrast η ranging from 1.5 to 4 to match with the typical amplitudes of the observed LPVs, and the exponents ($\alpha_0, \alpha_1, \alpha_2$) = (1, -14, 5).

The modelled LPVs are in good agreement with the observed ones, especially in Part IV where the best phase coverage was achieved in the observations, allowing for the clear detection of the double 'S' pattern due to the two CIR arms. In the modelled CIR patterns corresponding to Part I, we notice that we do not reproduce the apparent bump happening around the line centre around phase ~0.8. This could be due to the fact that the LI algorithm might have missed detecting a faint/small spot that would have caused a CIR arm responsible for that bump during Part I. However, it has to be noted and emphasized that this part of the observing campaign has the worst phase coverage, which plays an important role in the clear detection and characterization of the LPVs.

### 4.1.3 Variability in other lines

In view of the behaviour in the phased diagrams for He II λ4686, we explored the variability in other line profiles, some photospheric (e.g. He II λ5411) and some forming in the wind (Hα). Our CLEAN Fourier analyses across the lines detected the rotation period in most of them (Fig. 17), and even the first harmonic of the rotation period is detected in the Balmer lines. However, the phased dynamic difference diagrams for these lines, either by considering the entire time series or by splitting it into different parts, does not show any obvious pattern, which might be an indication on the weakness of the signal in these lines, requiring more S/N for their characterization in phased dynamic difference plots. Particularly, in the cases of photospheric lines, we do not detect any sign of the ~8.5 h periodicity reported by Baade (1986) and Reid & Howarth (1996), except a possible marginal detection in the blended He I λ6678+He II λ6683 line for which the noise level in its CLEAN periodogram is much higher than those of the other lines.

To further inspect the correlation between the variations observed by BRITE and the LPVs, we extracted the time series of the equivalent widths (EW) of He II λ4686, Hα and Hβ (Fig. 18). Not only is the behaviour of the phased EW variations and the phased light variations strongly correlated, but more importantly a noticeable phase delay is observed in the EW variations in the lines most sensitive to the wind (He II λ4686 and Hα) with respect to the photospheric light variations probed by BRITE. Table 6 lists the outcome of a cross-correlation analysis of the phased EW measurements with the phased BRITE light curve averaged in 0.04 phase bins. From these cross-correlation analyses, we conclude that, besides the noticeable phase delays in He II λ4686 and Hα with respect to the light variations probed by BRITE, there appears to be a real phase delay even between He II λ4686 and Hα, and between Hβ and the BRITE observations. More importantly, the phase delays observed in Fig. 18 and quantified by the cross-correlation analyses are indicative of the azimuthal difference $\Delta\varphi = \varphi_{CIR} - \varphi_s$ involved in equation (29), such that we can solve for r in that relation and thus extract the typical radial distance in the wind at which we probe the CIR arms through Hβ, Hα and He II λ4686. The resulting values of r are also listed in Table 6, in which the error bars are the consequences of accounting for the uncertainties on R, $v_\infty$, Ω and $\Delta\varphi$. It is worth noting that the values that
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Normalized CLEAN power spectra of the time series of difference profiles of some optical lines in ζ Pup, integrated over ±500 km s⁻¹ except for He II λ6560+Hα for which integration was performed over ±1000 km s⁻¹. The 1.78-d period due to the stellar rotation is detected in all these lines. The vertical red line points to the location of the first harmonic of the rotation period (here detected in the Balmer lines), while the vertical dotted green line indicates the location of the 5.075 ± 0.003-d period (no detection here), and the vertical black dashed line indicates the location of the 8.54 ± 0.054 h period (only a possible marginal detection in the noisy and blended He I λ6678+He II λ6683 lines).</figcaption></figure>

we find for *r* are in good agreement with the extent of the regions of line formation (up to ~2.5–3.5*R* for He II λ4686 and Hα as previously mentioned, and closer to the stellar surface for Hβ).

Thus, from all these considerations, we conclude that the 1.78-d cyclic intrinsic light variations observed by *BRITE* in ζ Pup, which we found to occur at the photosphere (Section 3.2) and are due to rotational modulation related to the presence of bright photospheric inhomogeneities (Section 3.4), are the photospheric drivers of the large-scale CIR arms in the stellar wind as traced by the cyclic variations that we detect in He II λ4686 (and Hα). This implies that our findings constitute the first observational evidence for a direct link between CIRs in the wind of an O-type star and their photospheric drivers.

<figure><img src="image_20.png" /><figcaption>**Figure 18.** First four panels from bottom to top: light variations of ζ Pup as observed by *BRITE* phased with the stellar rotation period, along with the phased time series of EWs of the Hβ, Hα and He II λ4686 line profiles of the star. The vertical scale on the panel for Hβ is exactly twice larger than the vertical scales on He II λ4686 and Hα. Upper panel: cross-correlation sequences of the EW variations with respect to the light variations observed by *BRITE*. The measured phase delays extracted through Gaussian fits to the first main cross-correlation peaks are indicated by the vertical markers (the Gaussian fits themselves are not overplotted for the sake of visibility).</figcaption></figure>

### 4.1.4 CIR/DAC recurrence time-scales

As mentioned in Section 1.1, the presence and blueward recurrent propagation of DACs observed in the absorption troughs of unsaturated resonance lines of most O-type stars is best interpreted as the spectroscopic manifestation of the presence of CIRs in their winds (Mullan 1984; Cranmer & Owocki 1996). Also, previous *IUE* monitoring of ζ Pup unveiled that the mean DAC recurrence period of the star was 19.23(45) h, so that when it was believed that the stellar rotation period is 5.1 d, the star would exhibit on
<table><caption>**Table 6.** Phase delays measured from the cross-correlations of the phased EW variations of the H β, H α and He II λ4686 lines of ζ Pup with respect to its phased light variations observed by *BRITE*. The values of the phase delays Δφ were assessed through Gaussian fits to the first main cross-correlation peaks, then translated via equation (29) into the radial extent of the CIR arm (measured from the centre of the star) probed by each line profile.</caption><thead><tr><th>Line</th><th>Δφ</th><th>r (R)</th></tr></thead><tbody><tr><td>H β</td><td>0.04 ± 0.03</td><td>1.7 ± 0.7</td></tr><tr><td>H α</td><td>0.11 ± 0.01</td><td>2.9 ± 0.9</td></tr><tr><td>He II λ4686</td><td>0.14 ± 0.02</td><td>3.5 ± 1.1</td></tr></tbody></table>

average 5–6 DACs per rotation cycle (fig. 1 of Howarth et al. 1995), a behaviour that deviates from that found in most O-type stars for which an average of two DACs per rotation cycle are generally observed (Kaper et al. 1999). Now our results that the true rotation period of ζ Pup is much shorter, $P = 1.78063(25)$ d, changes that picture. It is worth noting that the consecutive bumps observed in the phased *BRITE* light curves and in the EW variations of the wind emission lines have a separation that varies between 0.4 and 0.5 in phase. For instance, a separation of $Δφ = 0.45$ in rotational phase would correspond to two ‘events’ separated by a time interval of $Δt = 19.2308(27)$ h, which is compatible with the DAC recurrence time-scale of 19.23(45) h reported by Howarth et al. (1995), but then in that configuration the next ‘event’ would have to happen after $Δφ' = 0.55$, i.e. $Δt' = 23.5043(33)$ h later, such that the Fourier analysis of the time series over several cycles would still give rise to a peak at the first harmonic (21.3676(57) h), not at the ~19.23 h time-scale. Additionally, it has to be considered that Prinja et al. (1992) identified a DAC recurrence time-scale of ~15 h (no error bar available) from the examination of ~2.2 d of intensive observations within a 5.5-d contiguous *IUE* monitoring of ζ Pup in 1989. Thus, all the information available to date appears to point out that the average inter-DAC ‘period’ of ζ Pup has not always been strictly the same over the past three decades. In Fig. 19, we plot the known ‘periodicities’ that might be associated with CIRs/DACs in ζ Pup. The change in ‘periodicity’ seems gradual over about two decades, increasing from ~15 to 21 h, although rapidly at first, then asymptotic later. In this plot, we have also marked the location of the $16.67 ± 0.81$ and $16.90 ± 0.48$ h periodicities reported by Berghoefer et al. (1996) from a simultaneous X-ray and H α monitoring of ζ Pup spanning 11 d in 1991, as they appear to follow the global trend in the plot, although at this point the physical origin of these periodicities is not well established and more importantly has never been formally proven to be related to CIRs/DACs. Also, for the epoch of the *Coriolis/SMEI* observations, we have explicitly marked the location of the first harmonic of the 1.78-d signal, but one has to keep in mind that only the fundamental frequency has a high statistical significance in these data sets.

The global trend in Fig. 19 is unlikely related to a slowing down (especially if non-uniform in time) of the stellar rotation period: the mass-loss rate is simply too low to allow this and there is no plausible or conceivable internal mechanism to cause such a large change in rotation over such a short time. On the other hand, looking at the behaviour of the Sun in terms of rotation and spot activity, we see clear surface differential rotation from a period of ~25 d at the equator to ~35 d at the poles, with spot migration from mid to low latitude over an 11-yr cycle. Among late-type stars in

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Known ‘periodicities’ that might be related to CIRs/DACs in ζ Pup. Note: The ~15 h time-scale reported by Prinja et al. (1992) does not have an error bar. The results of the simultaneous ground-based photometry and H α spectroscopy in support of the *IUE* MEGA campaign observations remain unpublished, although they reveal a 19.7(1.1) h periodicity that is consistent with the results of the *IUE* observations.</figcaption></figure>

general, one does see both solar-type surface differential rotation (faster rotation at the equator: e.g. Donati & Collier Cameron 1997; Kovári et al. 2004; Petit et al. 2004; Roettenbacher et al. 2011) and anti-solar type of differential rotation (slower rotation at the equator: e.g. Hackman et al. 2001; Strassmeier, Kratzwald & Weber 2003; Harutyunyan et al. 2016). Inspired by this, we speculate that ζ Pup with its very rapid overall rotation could also exhibit differential rotation with the dominant spots wandering quickly first from higher latitudes with faster rotation to spots wandering slowly to lower latitudes with slower rotation at the bulging equator. In that scenario, the dominant spots observed at a given epoch of observation would be located around the same latitude, forming an ‘active latitude’, analogous to the Sun’s active latitudes that change over time to form the so-called Butterfly diagram of the Sun’s activity cycle. That being said, the surface maps that we obtained through the LI of the *Coriolis/SMEI* and *BRITE* light curves of ζ Pup also reveal the presence of spots at mid-latitudes. However, one important point that we emphasize here is that the LI algorithm as it stands currently allows for rigid rotation only, although it still can reveal differential rotation from the relative drifts in spot longitudes between different data sets (Roettenbacher et al. 2011). Nevertheless, on this note we leave further speculation in anticipation of future long-term monitoring of this ever-more interesting massive star.

## 4.2 Small-scale wind structures

The densest temporal coverage for the multisite ground-based spectroscopic monitoring of ζ Pup that we conducted in parallel with the *BRITE* observing run was achieved during ~5 d in 2015 February (HJD 2, 457, 059.94 – 2, 457, 065.10) when we had optimal longitude coverage from facilities at CTIO-SMARTS 1.5 m/CHIRON, SAAO 1.9 m/GIRAFFE and SASER. We took this opportunity to look for signatures of small-scale wind structures in He II λ4686, as previously found by Eversberg et al. (1998).
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Kinematics of outmoving wind clumps in ζ Pup as traced in He II λ4686, along with their photospheric drivers as observed by BRITE during the night of 2015 February 10/11. Top: average profile of the He II λ4686 line for the night. Middle: dynamic difference spectra from the average spectrum after removal of the effects of the 1.78-d cyclic variations due to large-scale wind structures (Section 4.2.1). Small horizontal ticks on the right-hand side of the diagram indicate the time sampling, with observations performed from Shenton Park Observatory (black), SAAO 1.9 m/GIRAFFE (green) and CTIO-SMARTS 1.5 m/CHIRON (red). Bottom: average scalogram from the Ricker wavelet-based CWTs of the spectra (in the spatial/wavelength domain). Left: light variations of ζ Pup recorded by BRITE in the blue (blue points) and red (red points) filters during the night of 2015 February 10/11. The light curves illustrated here are those free of the effects of the 1.78-d signal associated with rotational modulation (Section 3.5, Figs 11 and 12). Green squares represent the integrated intensity variations in the difference spectra within the velocity range ±1000 km s⁻¹, with the scaling factor C = −3 (equation 30, Section 4.2.3). Far left: portion of the scalogram of the residual BRITE light curve restricted to that specific night (taken from Fig. 12), showing that the dominant time-scale during that night was ~14–15 h. There seems to be a slight pseudo-period drift during the night, as already noted in Section 3.5.2. Right: strong correlation between the integrated intensity variations in He II λ4686 and the light variations observed by BRITE.</figcaption></figure>

### 4.2.1 Dynamic difference spectra

Our spectra have a typical S/N ~ 500 in the continuum close to He II λ4686, which is about a factor two lower than the S/N in the CFHT 3.6 m/Reticon spectra of ζ Pup in which Eversberg et al. (1998) found manifestations of wind clumps. Still, we detect signatures of wind clumps in our observations, as clearly revealed by the dynamic difference spectra in the main panel of Fig. 20 for the night of 2015 February 10/11 (HJD 2, 457, 064.03 – 065.10) during which the signatures are the most discernible: features in emission excess that appear at random both in Doppler velocity and in time, and propagate away from the line centre towards the red/blue edges while smearing out and dropping in intensity with time until they reach Doppler velocities ~±500 km s⁻¹ from where they seem to fade away. That the features seem to disappear beyond velocities ~±500 km s⁻¹ does not necessarily mean that the actual clumps causing these signatures physically vanish, as the LER in which we see their manifestations only spans a limited Doppler velocity range

that clearly does not reach the terminal wind speed in the case of He II λ4686 for ζ Pup which has $v_∞ = 2300 ± 100$ km s⁻¹. Also, another striking property that we notice in the dynamic difference spectra on the main panel of Fig. 20 is that the velocity widths are narrower for features located close to the line centre (e.g. those that appear at time ~5519.50) compared to those features that appear at the wings of the emission line, a behaviour also seen in other hot stars (Moffat 1994; Moffat et al. 1994; Lépine & Moffat 1999) and best interpreted as the manifestation of anisotropic velocity dispersion seen in Doppler projection.

We note that, at first, from a preliminary look at the dynamic difference spectra we wondered how the cyclic 1.78-d signal associated with large-scale wind structures influences the shape of the features related to clumps. The major difficulty here is the evolutive nature of the CIRs in time (ultimately as a consequence of the evolution of their photospheric drivers as discussed in Sections 3.4 and 4.1), but more importantly the ~5-d intensive observations
including the night of 2015 February 10/11, turn out to be right in the middle of the transition phase in Part III where we only see mild signature of CIRs in the phased diagrams, unlike in Part IV where the CIR signatures are more clearly defined. Nevertheless, we proceeded with the removal of the observed mild CIR patterns in Part III by adopting the same approach as we used for the removal of the 1.78-d signal in the light curves (Section 3.5), but this time performed for each Doppler velocity in the range ±1275 km s⁻¹: we phase-folded the time series spectra in Part III with the rotation period, rebinned them in phase, interpolated over the original phase sampling, subtracted the resulting smoothed pattern from the original phased time series, and then went back in the time domain to compare the resulting corrected time series with the original time series in order to check for possible overcorrections in which case we restarted again the whole procedure with another phase bin size. In general, during the comparison step we noted that the effect of the 1.78-d variability on the ~1.08-d time span of the 2015 February 10/11 night was a slight increase of intensity between times ~5519.70 and 5520.10. The phase bin size is the key parameter for controlling overcorrections and undercorrections, since too small phase bin sizes result in templates with relatively large noisy fluctuations (leading to overcorrection), while too large phase bin sizes yield too smooth templates (yielding undercorrection). We found that a phase bin size of 0.05 was a good compromise, yielding the dynamic difference spectra illustrated on the main panel of Fig. 20.

### 4.2.2 Wavelet analysis

Given the stochastic nature of the appearance of the spectroscopic features associated with the wind clumps, and also their morphology as small features in emission excess superimposed on the wind emission line, the wavelet formalism (here in the spatial/wavelength domain) is very appropriate for the extraction of the properties of the observed clumps, as was already performed in other investigations on spectroscopic observations of WR stars (Lepine 1994; Moffat 1994; Moffat et al. 1994; Lepine, Moffat & Henriksen 1996; Lépine & Moffat 1999). Also, unlike the case of the residual light curve on which we performed a wavelet analysis to characterize the time-scales of the features composing the observed stochastic signal intrinsic to ζ Pup (Section 3.5.2), our spectra have a constant dispersion in wavelength<fnref n="3" /> and thus are intrinsically equally sampled, such that the continuous formalism can be applied and the CWTs readily calculated. Still, the choice of an appropriate wavelet family to be used needs to be considered. Since the subpeaks that we want to characterize in the difference spectra are excesses in emission often accompanied by less strong side depressions (unlike the features in the light curve that have roughly equal bumps and dips between ±10 mmag), the simplest and most appropriate family of wavelets for this situation is the family of Ricker wavelets, also often called Mexican hat wavelets. The Ricker wavelet function is constructed from a normalization of the second derivative of a Gaussian function. It was also the type of family of wavelets adopted by the studies on observations of clumping in hot stars mentioned previously. Under these considerations, we calculated the CWTs

of the 14 spectra during the night of 2015 February 10/11 and extracted the mean scalogram which is depicted in the lower panel of Fig. 20, confirming that the substructures are mostly located between ~±500 km s⁻¹, but more importantly that the typical widths of the detected features as indicated in the scalogram are of the order of 10².⁰⁻¹⁰²¹ km s⁻¹ which is slightly higher but remains consistent with the value of ~10².⁰ km s⁻¹ found by Lépine & Moffat (1999) for the clumps in WR stars.

Also in order to fully compare the characteristics of the clumps of ζ Pup observed in our campaign with the clumps detected in the star at the epoch of the observations performed by Eversberg et al. (1998), we re-extracted the archival CFHT 3.6 m/Reticon optical spectra of the star taken during the nights of 1995 December 10/11 (HJD 2, 450, 061.94 – 062.15) and 1995 December 12/13 (HJD 2, 450, 063.95 – 064.17), and calculated the corresponding average Ricker wavelet-based scalograms, which are displayed in Fig. 21 below the corresponding dynamic difference spectra. The average scalograms for the CFHT 3.6 m/Reticon spectra indicate that the typical widths of the features are also ~10².⁰ km s⁻¹, with some isolated minor power detected at larger scales ~10².<fnref n="3" /> km s⁻¹. The presence of such broader features could be due to either nesting effects or contamination from large-scale wind structures, the latter being impossible to remove in the CFHT 3.6 m/Reticon spectra since the gap between the two nights of observations is of the order of one rotational cycle and each night only spanned ~5 h.

### 4.2.3 Correlated light variations and LPVs

An inspection of the behaviour of the stochastic light variations of ζ Pup as recorded by BRITE (Section 3.5) during the night of 2015 February 10/11 shows a relatively wide bump that reaches a maximum at time ~5519.70 (left-hand panel of Fig. 20), and the Morlet wavelet-based scalogram of the BRITE light curve indicates that the dominant pseudo-period during that night is ~14–15 h with a slight pseudo-period drift towards shorter time-scales during the night (far left-hand panel of Fig. 20). As already noted in Section 3.5.2, that pseudo-period drift, of unexplained origin, is part of a drift that starts from time ~5517.5 at time-scales ~16 h and finishes at time ~5521.0 at time-scales ~11 h. Nevertheless, in order to quantify the correlation between the stochastic light variations observed by BRITE in ζ Pup and the observed clumps in the wind, we assessed the integrated residual intensity variations in each difference spectrum over Doppler velocities ±1000 km s⁻¹:

$$ \mathcal{I}_k = \frac{C}{v_2 - v_1} \int_{v_1}^{v_2} I_r(v, t_k) dv \quad (30) $$

in which $I_r(v, t_k)$ is the residual intensity at Doppler velocity $v$ in the difference spectrum observed at time $t_k$, and the integration is performed over velocities ranging from $v_1 = -1000 \text{ km s}^{-1}$ to $v_2 = +1000 \text{ km s}^{-1}$. The constant $C$ is a negative normalization factor accounting for the comparison with the BRITE light variations which are expressed in magnitudes. The individual values of $\mathcal{I}_k$ for the 14 spectra in the night are displayed as green squares on the left-hand panel of Fig. 20 along with the residual BRITE light curves, taking $C = -3$. In view of that panel it appears that correlation exists between the two, an impression that is confirmed by the resulting Pearson correlation coefficient value $\rho = 0.8$ (right-hand panel of Fig. 20).

Since BRITE probes the light variations of ζ Pup at the level of its photosphere (Section 3.2), we conclude that the stochastic light variations detected by BRITE in the star, which we characterized in

<footnote n="3">³ In general, raw echelle spectra do not have a constant dispersion through the different orders, but are resampled with a constant dispersion during the process of merging of the orders and extraction of the one-dimensional wavelength-calibrated spectra.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Outmoving wind clumps in ζ Pup as first reported by Eversberg et al. (1998) from two nights of CFHT 3.6 m/Reticon observations (archival data) in 1995 December 10/11 (left) and 1995 December 12/13 (right). For a given night, the top panel illustrates the mean spectrum for the night, the middle panel shows the dynamic difference spectra from the average spectrum and the bottom panel depicts the average scalogram. The scalogram signatures of the clumps are the same as observed during the BRITE 2014–2015 observing run, with a dominant scale at a ~ 10<sup>2.0</sup>–10<sup>2.1</sup> km s<sup>-1</sup>.</figcaption></figure>

Section 3.5, are the manifestations of the photospheric drivers of the clumps in the stellar wind as traced in He II λ4686.

## 5 DISCUSSION AND CONCLUSION

The relatively short rotation period of ζ Pup, P = 1.78063(25) d, that we inferred from the photometric observations has a strong impact on our current understanding of the properties of the star. Most importantly, the now better constrained relatively low stellar inclination angle i ≃ 24° and high equatorial velocity (Section 3.4.1) completely change the previous beliefs that the star has a rotation period of ~5.1 d, is seen almost equator on and spins at $v_e \simeq v_e \sin i = 219 \pm 18$ km s<sup>-1</sup>. Moreover, as mentioned earlier (Section 3.4.1), the short rotation period of ζ Pup is in good agreement with the suspected rotational evolution of the star, with two proposed scenarios so far: either spin-up during a Roche lobe overflow (RLOF) phase in a WR+O close massive binary followed by ejection from the supernova explosion of the primary component leading to the current status of ζ Pup as a single fast rotator runaway (Vanbeveren et al. 1998), or a spin-up resulting from the merger of at least two massive stars during dynamical interactions within the Vela R2 R-association (Vanbeveren et al. 2009; Vanbeveren 2012; Pauldrach, Vanbeveren & Hoffmann 2012). In the former case, Vanbeveren et al. (1998) estimated that, with initial binary masses of 40 and 38 M<sub>☉</sub> and an initial orbital period of 4 d, the resulting orbital period of the system at the end of the RLOF phase is 2.6 d, in which case tidal effects would be very large and synchronize the components of the system such that the rotation period of the secondary would also be tuned to 2.6 d (fig. 14.15 in Vanbeveren et al. 1998, p. 232). This is of the same order of magnitude as the actual value of the rotation period that we inferred from the photometric observations in this investigations, P = 1.78063(25) d. The exact values of the initial binary masses implied by this post-RLOF rotation period depend on the exact values of the physical parameters involved during the RLOF, e.g. the exact amount of mass and

angular momentum loss from the binary. However, in that scenario, the runaway velocity acquired by the secondary after the supernova explosion (asymmetrical) of the primary has to satisfy (Vanbeveren et al. 1998, equations 11.4 and 11.17):

$$ \| \vec{v}_{\text{rw}} \| = \frac{M_{\text{WR}}}{\sqrt{M_{\text{WR}} + M_{\text{O}}}} \left(\frac{G}{A}\right)^{\frac{1}{2}} \quad (31) $$

in which the masses $M_{WR}$ and $M_O$, as well as the semi-major axis A are values just before the supernova explosion (here G denotes the gravitational constant). Assuming pre-supernova masses $M_{WR} = 10 M_{☉}$ and $M_O = 56 M_{☉}$, the runaway velocity corresponding to an orbital period of ~1.78 d would be ~107 km s<sup>-1</sup>, which is too high compared to the observed runaway velocity of the star ($\|\vec{v}_{rw}\| = 60.0 \pm 16.9$ km s<sup>-1</sup>), while only an orbital period >10 d would yield a runaway velocity that is more compatible with the observed one. The latter case remains a plausible scenario since in that situation tidal effects are very small and the rotation period of the gainer (which was spun-up to ~1.78 d during the RLOF phase) remains unchanged. Additionally, in order to satisfy the observed current mass of ζ Pup, $M = 56.1^{+14.5}_{-11.6} M_{☉}$ (despite the huge uncertainty on that value), the initial binary masses have to be much larger than those assumed in fig. 14.15 of Vanbeveren et al. (1998), e.g. 55 and 50 M<sub>☉</sub>. In that case, it has to be assumed that the 55 M<sub>☉</sub> primary ends its life with a supernova explosion that disrupted the binary, which is not entirely excluded as there is some evidence that black hole formation (at least for the lowest mass black holes) can be accompanied by a supernova explosion (e.g. Israelian et al. 1999). Under these circumstances, the updated scenario is illustrated in Fig. 22. We emphasize that we performed these calculations only to update this scenario to satisfy the most pertinent values of the stellar parameters of ζ Pup, but at this point no strong constraints is available to help decide which of the two proposed evolutionary scenarios is the most plausible one.

With respect to the physical mechanism at the origin of the photospheric bright spots that give rise to the rotationally modulated light
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Evolutionary history of $\zeta$ Pup in the scenario in which the star could have been the secondary component of a massive binary and would have been spun-up during a RLOF phase of the primary component (fourth row from the top). The latter later underwent a supernova explosion, leading to the current fast rotator runaway status of the $\zeta$ Pup. Figure adapted from Vanbeveren et al. (1998), updated with a rotation period of 1.78 d, and satisfying the observed runaway velocity $\|\vec{v}_{rw}\| = 60.0 \pm 16.9$ km s$^{-1}$ and the stellar mass $M = 56.1_{-11.6}^{+14.5} M_{\odot}$ listed in Table 1.</figcaption></figure>

variations observed by *BRITE* and *Coriolis/SMEI*, one possible candidate is the surface emergence of small-scale magnetic fields generated through dynamo action within a subsurface convection zone associated with the iron-opacity bump at $T \sim 150\,000$ K (FeCZ: Cantiello et al. 2009; Cantiello & Braithwaite 2011). The detection of such small-scale magnetic fields requires high-resolution, high-S/N circular spectropolarimetric observations imposed by the fast rotation of $\zeta$ Pup. All previous spectropolarimetric observations of the star ended up with no detection, with an upper limit of 121 G for a dipolar field strength, and longitudinal field error bar of 21 G (David-Uraz et al. 2014). Also since most of the spots that we detect during the *BRITE* and *Coriolis/SMEI* observing runs lie near the equator and are separated by $\sim 144^\circ - 180^\circ$ in longitude, their magnetic signature (if they are related to small-scale magnetic fields) might cancel out in low-resolution mode as was the case of the spectropolarimetric monitoring of $\zeta$ Pup performed by Hubrig et al. (2016) who, despite no magnetic detection and no period detection, found marginal signs of modulation of the longitudinal magnetic field measurements when they phase-folded their data with the 1.78-d period (rather than the 5.1-d period). More importantly, as mentioned in Section 1.1 previous photometric monitoring of O stars revealed light variations best interpreted as signs of rotational modulation possibly due to surface inhomogeneities, with only the case of $\xi$ Per showing potential link to its CIR/DAC behaviour (Ramiaramanantsoa et al. 2014), but our present findings that the bright localized photospheric spots observed by *BRITE* and *Coriolis/SMEI* in $\zeta$ Pup are the photospheric drivers of its CIRs/DACs constitutes the first observational evidence for a link between large-scale wind structures in an O-type star and their photospheric origin. This reinforces the need to observe more O-type stars with known CIR/DAC recurrence time-scale (e.g. those studied by Kaper et al. 1999) on a long-term through high-precision photometric monitoring and contemporaneous high-S/N mid/high-resolution spectroscopic/spectropolarimetric monitoring, because the near-ubiquity of CIRs/DACs in O-type stars implies that their photospheric drivers must also be a common phenomenon.

Lastly, the detection of the stochastically excited (but coherent for several hours) light variations intrinsic to $\zeta$ Pup during the *BRITE* and *Coriolis/SMEI* observations is as important as the cyclic variations related to rotational modulation itself, since we also find for the first time that those randomly excited (but coherent for several hours) light variations are the photospheric drivers of stellar wind clumps. Although the physical nature of these stochastic intrinsic variations at the photosphere remains unknown, there is a possibility that such variations arise from randomly excited stellar oscillations with finite lifetimes such as IGWs, generated within an FeCZ (Cantiello et al. 2009) or even from the convective core in the case of IGWs as found by Aerts & Rogers (2015). Nevertheless, the fact that these stochastic photospheric variations are linked to the variability associated with wind clumps suggests that the latter are formed even at the very base of the wind. Previous X-ray observations of $\zeta$ Pup suggested that only a large number of clumps in the wind ($N > 10^5$) are able to reproduce the observed X-ray variability of the star (Nazé et al. 2013). Since the variability scales with the inverse of the square root of the number of clumps involved, we can infer by the optical line flux variability that we observe ($\sim 6$ per cent) that it would arise from the contribution of $N \simeq 300$ clumps. The question therefore remains as whether those clumps are of equal size or follow a power-law distribution in which there is a relatively small number of strongly emitting clumps which would dominate the variability. The latter case fits naturally into the notion of the scaling laws of turbulence (compressible for hot winds) with a
power law of negative slope, making few large clumps along with a rapidly increasing number of smaller ones each of which emits less and less as the size is diminished (e.g. Moffat 1994). This is further supported by the fact that turbulence is a very common phenomenon in astrophysics, e.g. large-scale structures in the Universe, the inter-galactic medium (Ly-α clouds towards quasars), galaxies, cloudlets in giant molecular clouds, the discrete nature of stars themselves (as implied by the IMF) or surface granulation in the Sun itself. Also, the fairly clear and distinct observed spectral emission-line excess streaks over several hours must be spatially confined and not due to e.g. a sum of smaller random perturbations anywhere in the wind. However, it has also been shown that the first scenario, that is variability that arises from a large number of clumps of equal size, cannot be ruled out (e.g. Davies, Vink & Oudmaijer 2007). To further test the hypothesis of turbulence, high-S/N high-resolution spectroscopic observations of a WR star with strong emission line and high terminal wind speed (e.g. up to ~5000 km s⁻¹) needs to be performed, the best candidates being one of the WO2 stars in the Milky Way. But given the faintness of these objects, observations with the upcoming generation of extremely large telescopes such as the Thirty Metre Telescope would be the most ideal.

## ACKNOWLEDGEMENTS

We gratefully acknowledge fruitful discussions with D. John Hillier, Yaël Nazé and Radosław Smolec in relation to this investigation. We also thank the reviewer, Richard Townsend, for his insights and for providing us with useful constructive remarks. TR acknowledges support from the Canadian Space Agency (CSA) grant FAST. AFJM is grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Fonds de Recherche du Québec – Nature et Technologies for general financial aid, and to CSA for funding the training of highly qualified personnel associated with the BRITE project. NDR acknowledges postdoctoral support by the University of Toledo and by the Helen Luedtke Brooks Endowed Professorship. API acknowledges support from the National Science Centre grant no. 2016/21/B/ST9/01126. APo acknowledges National Science Centre (NCN) grant 2016/21/D/ST9/00656. GH acknowledges support from NCN grant 2015/18/A/ST9/00578. The Polish contribution to the BRITE mission is supported by the Polish NCN grant 2011/01/M/ST9/05914. GAW and SMR acknowledge Discovery Grant support from NSERC. KZ acknowledges support from the Austrian Fonds zur Förderung derwissenschaftlichen Forschung (project V431-NBL). This investigation is based on data collected by the BRITE-Constellation satellite mission, designed, built, launched, operated and supported by the Austrian Research Promotion Agency (FFG), the University of Vienna, the Techni-cal University of Graz, the CSA, the University of Toronto Institute for Aerospace Studies, the Foundation for Polish Science and Technology and NCN. This project also uses observations made from the SAAO. The professional authors of this paper are grateful to the amateur astronomers of the SASER team (also co-authors of this paper), who invested personal time and enthusiasm in this project. Part of the research leading to these results has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 670519: MAMSIE).

## REFERENCES

* Abt A., 1952, AJ, 57, 158

* Aerts C., Rogers T. M., 2015, ApJ, 806, L33

* Aerts C. et al., 2013, A&A, 557, A114

* Aerts C. et al., 2017, A&A, 602, A32

* Auvergne M. et al., 2009, A&A, 506, 411

* Baade D., 1986, NATO Adv. Sci. Inst. Ser. C, 169, 465

* Baglin A. et al., 2006, in Fridlund M., Baglin A., Lochard J., Conroy L., eds, ESA SP-1306, The CoRoT Mission Pre-Launch Status – Stellar Seismology and Planet Finding. ESA, Noordwijk, p. 33

* Balona L. A., 1992, MNRAS, 254, 403

* Balona L. A., 2016, MNRAS, 457, 3724

* Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339

* Beals C. S., 1953, Publ. Dom. Astrophys. Obs. Victoria, 9, 1

* Béky B., Holman M. J., Kipping D. M., Noyes R. W., 2014, ApJ, 788, 1

* Berghoefer T. W., Baade D., Schmitt J. H. M. M., Kudritzki R.-P., Puls J., Hillier D. J., Pauldrach A. W. A., 1996, A&A, 306, 899

* Bernhard K., Hümmerich S., Otero S., Paunzen E., 2015a, A&A, 581, A138

* Bernhard K., Hümmerich S., Paunzen E., 2015b, Astron. Nachr., 336, 981

* Blomme R., van de Steene G. C., Prinja R. K., Runacres M. C., Clark J. S., 2003, A&A, 408, 715

* Blomme R. et al., 2011, A&A, 533, A4

* Borucki W. J. et al., 2010, Science, 327, 977

* Bouret J.-C., Hillier D. J., Lanz T., Fullerton A. W., 2012, A&A, 544, A67

* Breger M. et al., 1993, A&A, 271, 482

* Breger M. et al., 1999, A&A, 349, 225

* Briquet M. et al., 2011, A&A, 527, A112

* Bruntt H., 2007, Commun. Asteroseismol., 150, 326

* Bruntt H. et al., 2009, MNRAS, 396, 1189

* Buysschaert B. et al., 2015, MNRAS, 453, 89

* Buysschaert B. et al., 2017, A&A, 602, A91

* Buzasi D. L., 2000, in Teixeira T. C., Bedding T. R., eds, Experiment Design and Data Reduction for Seismology: Lessons Learned from WIRE. Aarhus Universitet, Aarhus, Denmark, p. 9

* Buzasi D. et al., 2000, ApJ, 532, L133

* Cantiello M., Braithwaite J., 2011, A&A, 534, A140

* Cantiello M. et al., 2009, A&A, 499, 279

* Clarke D., 2003, A&A, 407, 1029

* Craig I. J. D., Brown J. C., 1986, Research supported by SERC. Bristol, England and Boston, MA, Adam Hilger, Ltd., 159 p.

* Cranmer S. R., Owocki S. P., 1996, ApJ, 462, 469

* Davenport J. R. A., Hebb L., Hawley S. L., 2015, ApJ, 806, 212

* David-Uraz A. et al., 2014, MNRAS, 444, 429

* Davies B., Vink J. S., Oudmaijer R. D., 2007, A&A, 469, 1045

* Davis J., Morton D. C., Allen L. R., Hanbury Brown R., 1970, MNRAS, 150, 45

* De Cat P., Aerts C., 2002, A&A, 393, 965

* Degroote P. et al., 2010, A&A, 519, A38

* Degroote P. et al., 2011, A&A, 536, A82

* Dessart L., Chesneau O., 2002, A&A, 395, 209

* Dessart L., Owocki S. P., 2003, A&A, 406, L1

* Donati J.-F., Collier Cameron A., 1997, MNRAS, 291, 1

* Eversberg T., Lépine S., Moffat A. F. J., 1998, ApJ, 494, 799

* Eyles C. J. et al., 2003, Sol. Phys., 217, 319

* Feldmeier A., 1995, A&A, 299, 523

* Ferraz-Mello S., 1981, AJ, 86, 619

* Gies D. R., Williams S. J., Matson R. A., Guo Z., Thomas S. M., Orosz J. A., Peters G. J., 2012, AJ, 143, 137

* Gontcharov G. A., 2006, Astron. Lett., 32, 759

* Hackman T., Jetsu L., Tuominen I., 2001, A&A, 374, 171

* Hamann W.-R., Gräfener G., 2003, A&A, 410, 993

* Hamann W.-R., Gräfener G., 2004, A&A, 427, 697

* Harmon R. O., Crews L. J., 2000, AJ, 120, 3274

* Harutyunyan G., Strassmeier K. G., Künstler A., Carroll T. A., Weber M., 2016, A&A, 592, A117

* Hervé A., Rauw G., Nazé Y., 2013, A&A, 551, A83

* Hillier D. J., Bouret J.-C., Lanz T., Busche J. R., 2012, MNRAS, 426, 1043

* Howarth I. D., 1992, in Drissen L., Leitherer C., Nota A. eds, ASP Conf. Ser. Vol. 22, Nonisotropic and Variable Outflows from Stars. Astron. Soc. Pac., San Francisco, p. 155

* Howarth I. D., Prinja R. K., 1989, ApJS, 69, 527

* Howarth I. D., Stevens I. R., 2014, MNRAS, 445, 2878
Howarth I. D., Prinja R. K., Massa D., 1995, ApJ, 452, L65

Howarth I. D., Siebert K. W., Hussain G. A. J., Prinja R. K., 1997, MNRAS, 284, 265

Howarth I. D., Goss K. J. F., Stevens I. R., Chaplin W. J., Elsworth Y., 2014, MNRAS, 440, 1674

Howell S. B. et al., 2014, PASP, 126, 398

Hubrig S., Kholtygin A., Ilyin I., Schöller M., Oskinova L. M., 2016, ApJ, 822, 104

Hundhausen A. J., 1972, in Sonett C. P., Coleman P. J., Wilcox J. M., eds, Solar Wind, Vol. 308. NASA Special Publication, Washington, DC, 393

Ignace R., Gayley K. G., 2003, MNRAS, 341, 179

Ignace R., Hendry M. A., 2000, ApJ, 537, L131

Ignace R., Hubrig S., Schöller M., 2009, AJ, 137, 3339

Ignace R., St-Louis N., Proulx-Giraldeau F., 2015, A&A, 575, A129

Israelian G., Rebolo R., Basri G., Casares J., Martín E. L., 1999, Nature, 401, 142

Kaper L., Henrichs H. F., Nichols J. S., Telting J. H., 1999, A&A, 344, 231

Kipping D. M., 2012, MNRAS, 427, 2487

Kolbin A. I., Shimansky V. V., 2014, Astrophys. Bull., 69, 179

Kolbin A. I., Shimansky V. V., Sakhibullin N. A., 2013, Astron. Rep., 57, 548

Kolbin A. I., Sakhibullin N. A., Gabdeev M. M., 2015, Adv. Space Res., 55, 808

Kovári Z., Strassmeier K. G., Granzer T., Weber M., Oláh K., Rice J. B., 2004, A&A, 417, 1047

Krtička J., Mikulášek Z., Zverko J., Žižkovský J., 2007, A&A, 470, 1089

Kuschnig R., Weiss W. W., Gruber R., Bely P. Y., Jenkner H., 1997, A&A, 328, 544

Lamers H. J. G. L. M., Cerruti-Sola M., Perinotto M., 1987, ApJ, 314, 726

Lamers H. J. G. L. M., Maeder A., Schmutz W., Cassinelli J. P., 1991, ApJ, 368, 538

Lenz P., Breger M., 2005, Commun. Asteroseismol., 146, 53

Lepine S., 1994, Ap&SS, 221, 371

Lépine S., Moffat A. F. J., 1999, ApJ, 514, 909

Lépine S., Moffat A. F. J., 2008, AJ, 136, 548

Lepine S., Moffat A. F. J., Henriksen R. N., 1996, ApJ, 466, 392

Lindegren L., Lammers U., Hobbs D., O'Mullane W., Bastian U., Hernández J., 2012, A&A, 538, A78

Loumos G. L., Deeming T. J., 1978, Ap&SS, 56, 285

Lucy L. B., Solomon P. M., 1970, ApJ, 159, 879

Lucy L. B., White R. L., 1980, ApJ, 241, 300

Lüftinger T. et al., 2010, A&A, 509, A43

Maíz-Apellániz J., Walborn N. R., Galué H. Á., Wei L. H., 2004, ApJS, 151, 103

Martín-Fleitas J. et al., 2014, Proc. SPIE, 9143, 91430Y

Martins F., Schaerer D., Hillier D. J., 2005, A&A, 436, 1049

Massa D. et al., 1995, ApJ, 452, L53

Mazzali P. A., 1990, A&A, 238, 191

McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203

Moffat A. F. J., 1994, Rev. Mod. Astron., 7, 51

Moffat A. F. J., Michaud G., 1981, ApJ, 251, 133

Moffat A. F. J., Drissen L., Lamontagne R., Robert C., 1988, ApJ, 334, 1038

Moffat A. F. J., Lepine S., Henriksen R. N., Robert C., 1994, Ap&SS, 216, 55

Moffat A. F. J. et al., 1998, A&A, 331, 949

Moskalik P., Buchler J. R., 1990, ApJ, 355, 590

Mullan D. J., 1984, ApJ, 283, 303

Natale G., Marconi M., Bono G., 2008, ApJ, 674, L93

Nazé Y., Oskinova L. M., Gosset E., 2013, ApJ, 763, 143

Oskinova L. M., Todt H., Ignace R., Brown J. C., Cassinelli J. P., Hamann W.-R., 2011, MNRAS, 416, 1456

Owocki S. P., Castor J. I., Rybicki G. B., 1988, ApJ, 335, 914

Pablo H. et al., 2015, ApJ, 809, 134

Pablo H. et al., 2016, PASP, 128, 125001

Pablo H. et al., 2017, MNRAS, 467, 2494

Pauldrach A. W. A., Vanbeveren D., Hoffmann T. L., 2012, A&A, 538, A75

Petit P. et al., 2004, MNRAS, 351, 826

Petrenz P., Puls J., 1996, A&A, 312, 195

Pigulski A. et al., 2016, A&A, 588, A55

Pollard K. R., Cottrell P. L., Kilmartin P. M., Gilmore A. C., 1996, MNRAS, 279, 949

Popowicz A., 2016, Proc. SPIE, 9904, 99041R

Popowicz A. et al., 2017, A&A, 605, A26

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 2007, Numerical Recipes, 3rd edn. Cambridge Univ. Press, Cambridge

Prinja R. K., 1988, MNRAS, 231, 21P

Prinja R. K. et al., 1992, ApJ, 390, 266

Puls J., Vink J. S., Najarro F., 2008, A&AR, 16, 209

Ramiaramanantsoa T. et al., 2014, MNRAS, 441, 910

Reeve D. C., Howarth I. D., 2016, MNRAS, 456, 1294

Reid A. H. N., Howarth I. D., 1996, A&A, 311, 616

Roberts D. H., Lehar J., Dreher J. W., 1987, AJ, 93, 968

Roettenbacher R. M., Harmon R. O., Vutisalchavakul N., Henry G. W., 2011, AJ, 141, 138

Roettenbacher R. M., Monnier J. D., Harmon R. O., Barclay T., Still M., 2013, ApJ, 767, 60

Roettenbacher R. M., Kane S. R., Monnier J. D., Harmon R. O., 2016, ApJ, 832, 207

Rybicki G. B., Hummer D. G., 1983, ApJ, 274, 380

Saio H., 2011, MNRAS, 412, 1814

Salpeter E. E., 1955, ApJ, 121, 161

Sander A., Shenar T., Hainich R., Gímenez-García A., Todt H., Hamann W.-R., 2015, A&A, 577, A13

Schumann J. D., Seggewiss W., 1975, in Proc. Symp. 67, Variable Stars and Stellar Evolution. D. Reidel Publishing Co., Dordrecht, p. 299

Shenar T. et al., 2015, ApJ, 809, 135

Shibahashi H., Osaki Y., 1981, PASJ, 33, 427

Sota A., Maíz Apellániz J., Morrell N. I., Barbá R. H., Walborn N. R., Gamen R. C., Arias J. I., Alfaro E. J., 2014, ApJS, 211, 10

Sprague T. B., 1990, PhD thesis, Western Michigan Univ.

St-Louis N., Tremblay P., Ignace R., 2017, preprint (arXiv:1710.09724)

Strassmeier K. G., Kratzwald L., Weber M., 2003, A&A, 408, 1103

Szabó R. et al., 2010, MNRAS, 409, 1244

Todt H., Sander A., Hainich R., Hamann W.-R., Quade M., Shenar T., 2015, A&A, 579, A75

Tody D., 1986, Proc. SPIE, 627, 733

Tody D., 1993, in Hanisch R. J., Brissenden R. J. V., Barnes J., eds., ASP Conf. Ser. Vol. 52, Astronomical Data Analysis Software and Systems II Astron. Soc. Pac., San Francisco, p. 173

Tokovinin A., Fischer D. A., Bonati M., Giguere M. J., Moore P., Schwab C., Spronck J. F. P., Szymkowiak A., 2013, PASP, 125, 1336

Twomey S., 1977, Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements. Elsevier, Amsterdam

Vanbeveren D., 2012, in Drissen L., Robert C., St-Louis N., Moffat A. F. J., eds, Proceedings of a Scientific Meeting in Honor of Anthony F. J. Moffat, ASP Conf. Ser., Vol. 465. Astron. Soc. Pac., San Francisco, p. 342

Vanbeveren D., van Rensbergen W., De Loore C., 1998, in Vanbeveren D., van Rensbergen W., De Loore C., eds, Astrophysics and Space Science Library, 232, The brightest binaries. Springer, The Netherlands

Vanbeveren D., Belkus H., van Bever J., Mennekens N., 2009, Ap&SS, 324, 271

Walker G. et al., 2003, PASP, 115, 1023

Walker G. A. H. et al., 2005, ApJ, 623, L145

Weiss W. W. et al., 2014, PASP, 126, 573

Weiss W. W. et al., 2016, A&A, 588, A54

Wild W. J., 1989, PASP, 101, 844

Wild W. J., 1991, ApJ, 368, 622

## APPENDIX A: EVOLUTION OF THE 1.78-D SIGNAL DURING THE *Coriolis/SMEI* OBSERVING RUN
<figure><img src="image_25.png" /><figcaption>**Figure A1.** Same as Fig. 7, but for the *Coriolis/SMEI* observing run in 2005–2006.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure A2.** Same as Fig. 7, but for the *Coriolis/SMEI* observing run in 2004–2005.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure A3.** Same as Fig. 7, but for the *Coriolis/SMEI* observing run in 2003–2004.</figcaption></figure>

## APPENDIX B: LIGHT CURVE INVERSION: SURFACE MAPS OF ζ PUP DURING THE *Coriolis/SMEI* SEASONAL OBSERVING RUNS

<figure><img src="image_28.png" /><figcaption>**Figure B1.** LI: mapping the photosphere of ζ Pup as observed by *Coriolis/SMEI* in 2005–2006, for $T^{(s)} = 42.5$ kK and a stellar inclination angle $i = 24^\circ$. Time increases upwards.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Figure B2.** LI: mapping the photosphere of ζ Pup as observed by Coriolis/SMEI in 2004–2005, for $T^{(s)}$ = 42.5 kK and a stellar inclination angle $i$ = 24°. Time increases upwards.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure B3.** LI: mapping the photosphere of ζ Pup as observed by Coriolis/SMEI in 2003–2004, for $T^{(s)} = 42.5$ kK and a stellar inclination angle $i = 24^\circ$. Time increases upwards.</figcaption></figure>

<footnote n="1">¹Département de Physique, Université de Montréal, CP 6128, Succursale Centre-Ville, Montréal, Québec H3C 3J7, Canada</footnote>

<footnote n="2">²Centre de Recherche en Astrophysique du Québec (CRAQ), Canada</footnote>

<footnote n="3">³Department of Physics and Astronomy, Ohio Wesleyan University, Delaware, OH 43015, USA</footnote>

<footnote n="4">⁴East Tennessee State University, Department of Physics and Astronomy, Johnson City, TN 37614, USA</footnote>

<footnote n="5">⁵Astrophysical Institute, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium</footnote>

<footnote n="6">⁶Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Str. 24/25, D-14476 Potsdam, Germany</footnote>

<footnote n="7">⁷Ritter Observatory, Department of Physics and Astronomy, The University of Toledo, Toledo, OH 43606-3390, USA</footnote>

<footnote n="8">⁸Department of Physics and Astronomy, University College London, Gower St, London WC1E 6BT, UK</footnote>

<footnote n="9">⁹School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="10">¹⁰Schnörringen Telescope Science Institute, Waldbröl, Germany</footnote>

<footnote n="11">¹¹Instytut Astronomiczny, Uniwersytet Wrocławski, Kopernika 11, PL-51-622 Wrocław, Poland</footnote>

<footnote n="12">¹²Instytut Automatyki, Politechnika Śląska, Akademicka 16, PL-44-100 Gliwice, Poland</footnote>

<footnote n="13">¹³Institut für Astrophysik, Universität Wien, Türkenschanzstrasse 17, A-1180 Wien, Austria</footnote>

<footnote n="14">¹⁴Nicolaus Copernicus Astronomical Center, Bartycka 18, PL-00-716 Warsaw, Poland</footnote>

<footnote n="15">¹⁵LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="16">¹⁶Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium</footnote>

<footnote n="17">¹⁷Department of Physics, Royal Military College of Canada, PO Box 17000, Station Forces, Kingston, Ontario K7K 7B4, Canada</footnote>

<footnote n="18">¹⁸Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario M5S 3H4, Canada</footnote>
This paper has been typeset from a TeX/LA TeX file prepared by the author.

<footnote n="19">¹⁹ Institut für Astro- und Teilchenphysik, Universität Innsbruck, Techniker-
strasse 25/8, A-6020 Innsbruck, Austria</footnote>

<footnote n="20">²⁰International Centre for Radio Astronomy Research, The University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="21">²¹SASER Team, Domain Observatory, 269 Domain Road, South Yarra, Vic 3141, Australia</footnote>

<footnote n="22">²²SASER Team, Dogsheaven Observatory, SMPW Q25 CJ1 LT10B, Brasilia, Brazil</footnote>

<footnote n="23">²³SASER Team, Latham Observatory, Australia</footnote>

<footnote n="24">²⁴SASER Team, R. F. Joyce Observatory, Canterbury Astronomical Society,
West Melton, New Zealand</footnote>

<footnote n="25">²⁵SASER Team, Mirranook Observatory, Booroolong Road, Armidale, NSW.
Australia</footnote>

<footnote n="26">²⁶Gemini Observatory, Northern Operations Center, 670 N. A'ohoku Place,
Hilo, HI 96720, USA</footnote>

<footnote n="27">²⁷South African Astronomical Observatory, PO Box 9, Observatory, 7935 Cape, South Africa</footnote>

<footnote n="28">²⁸Southern African Large Telescope, PO Box 9, Observatory, 7935 Cape, South Africa</footnote>

<footnote n="29">²⁹Eureka Scientific, Inc., 2452 Delmer Street, Oakland, CA 94602, USA</footnote>

<footnote n="30">³⁰Université de Toulouse, UPS-OMP, IRAP, Toulouse, France</footnote>

<footnote n="31">³¹CNRS, IRAP, 14, avenue Edouard Belin, F-31400 Toulouse, France</footnote>