A Study of the Stochastic Photometric Variability in the Winds of Galactic Wolf–Rayet Stars

The BRITE-Constellation (Weiss et al. 2014; Pablo et al. 2016) consists of five nanosatellites, each housing a 35 mm KAI-11002M CCD detector fed by a 30 mm diameter f/2.3 telescope through either a blue (390–460 nm) or a red (545–695 nm) filter. They are named BRITE-Austria (BAb), Uni-BRITE (UBr), BRITE-Heweliusz (BHr), BRITE-Lem (BLb) and BRITE-Toronto (BTr), the last letter of the abbreviations denoting the filter ("b" for blue and "r" for red). All the satellites were launched into low-Earth orbits with orbital periods of ∼100 minutes. With a ∼24° × 20° effective field of view (FOV) and a resolution of 27'' per pixel, each component of the BRITE-Constellation performs the simultaneous monitoring of 15 to 30 stars brighter than V ∼ 6. A given field is observed typically over a ∼6 month time period. As much as possible, at least two satellites equipped with different filters are set to monitor the field to provide dual-band observations.

After an initial extraction with the BRITE data reduction pipeline (Popowicz et al. 2017), the data are further inspected to ensure the removal of charge-transfer inefficiencies and hot pixels (Pablo et al. 2016). Then the data are made available in the BRITE archive. However, at this stage, the photometry may still be correlated with many satellite parameters, which are quantified by a variable in the data file depending on the data release version. The BRITE data sets used in this analysis are all from the fifth data release (DR5), and the intrinsic parameters against which the data can be decorrelated are presented in Table A.1 of Pablo et al. (2016). The decorrelation procedure we applied to our observations is the same as that presented in Appendix A of Pigulski et al. (2016).

Since the presence of a bright companion can significantly affect the observed variability, we did not include any known binary system in our sample except for the observations of the binary system WR22, since Lenoir-Craig et al. (2021) determined, through the analysis of photometric observations of the eclipse shape, the O-star component to be of spectral class O9V and thus of contributing only ∼8% of the WR flux to the luminosity of the system according to the stellar parameter calibrations of Martins et al. (2005). Seventeen observing runs including WR stars in the targeted fields were performed between 2014 and 2019 by BRITE. Of the WR stars observed, only WR6, WR24, and WR40 are single stars. We then checked the Gaia DR2 archive to ensure the absence of any bright visual companion within a 10 pixel radius of these stars. In Table 1 we list the name of the star and the year in which it was observed (if this stretched over two years, we refer to only the longest part), its spectral type from the online WR catalog maintained by P. Crowther, ¹⁰ the time interval over which it was observed, the satellite used, and the corresponding field identification. An example of a BRITE lightcurve is shown in Figure 1, where we show mean fluxes per orbit measurements for the WN6ha star WR24 as a function of the Heliocentric Julian Date. The overall median flux value has been subtracted from all data values.

Zoom In Zoom Out Reset image size

The main advantage of BRITE photometry is the generally relatively long duration of the observations of each field, which ranges from 62 to 165 days for the data used in the present study. In particular, the longer among these runs allow for a much more accurate probing of the variability in the low-frequency range.

The TESS has a 90° × 24° FOV with its four main cameras and performs observations of stars in the 600–1000 nm range (Ricker et al. 2015). Its average observing interval on a specific field is ∼27 days, with either a 2 or 30 minute cadence for data collected during its prime mission. A 20 s cadence was added in the extended mission that started on 2020 July 5, along with a 10 minute cadence replacing the previous 30 minute cadence. Its accurate photometry makes this satellite ideal for observing short-term variability in the winds of WR stars in the ν ∈ [0.03, 350] days⁻¹ domain. TESS data are publicly available on the Mikulski Archive for Space Telescopes (MAST) and can be accessed through various means. During the refereeing process of this paper, we became aware of another publication (Nazé et al. 2021) using a similar set of TESS observations to the one we are using in this study. They also report on the presence of red noise and discuss its characteristics.

We have retrieved all TESS data sets of single WR stars with a Gaia magnitude over 12 that were available on 2021 August 10. The details of the retrieved data sets are presented in Table 2, where we list the stars' names, their spectral types from the online WR catalog of P. Crowther (see references therein), the TESS sector in which the targets were located, the time intervals over which they were observed, the cadences used, and the fitted parameters of the semi-Lorentzian model described in Section 3.2. For each data set, we first checked for bright neighbors that could have contributed to the observed flux with the Gaia DR2 archive, and then cross-matched these with the Target Pixel Files (TPFs) also obtained through the Lightkurve package from the MAST, which shows how many pixels with length 21'' the source and neighbors occupy in the field. Stars with a bright, close-by (i.e., within a radius of less than 63'' from the target) visual neighbor were removed from our TESS sample, as were stars with bright visual neighbors further away that had a significant part of their flux smeared inside the 63'' radius around the target. We also did not include WR31b in our sample, as Groh et al. (2009) observed its spectral type varying from WN11h at minimum apparent magnitude of its S-Dor cycles to a type resembling an extreme A-type hypergiant at maximum visual magnitude, thus making it a different object than the cWR and WNh stars in our sample.

The 2 minute TESS lightcurves, which have been detrended by the TESS data pipeline following Jenkins et al. (2016), were obtained from the MAST using the Lightkurve package for Python (Lightkurve Collaboration et al. 2018), while the 10 and 30 minute TESS lightcurves were extracted from the TPFs. Aperture photometry with a custom mask for the source was performed on each data set, with a flux threshold usually in the range of 7–10 median absolute deviations (MADs) above the median flux, but sometimes lower for fainter sources or higher in the case where another bright star is present in the field. A background mask was defined by pixels with a median flux value of 10⁻⁴ MADs or less, and the background signal was subsequently subtracted from the source. A principal component analysis was then performed to remove the remaining instrumental noise and systematics from each data set, using the Lightkurve Regression Corrector method. An example of the resulting photometry (with the median value subtracted) is shown as a function of the Barycentric TESS Julian Date (BTJD) in Figure 2, where we present the 30 minute cadence TESS lightcurve of the WN8h star WR40.

Zoom In Zoom Out Reset image size

The MOST satellite (Matthews 2003, 2004; Walker et al. 2003) housed a 15 cm optical telescope feeding a CCD photometer through a single custom broadband (400–700 nm) filter. MOST followed a polar Sun-synchronous orbit at 820 km altitude, enabling it to monitor stars in its equatorial continuous viewing zone for up to 8 weeks without interruption with a cadence of ∼100 minutes.

Our sample includes six single WR stars observed with MOST, which are listed in Table 3 along with their spectral types, observing interval, and fitted parameters of the semi-Lorentzian model described in Section 3.2. However, only WR92, 103, 115, and 120 still had data files available with parameters related to the satellite's onboard instruments, allowing us to carry out detrending procedures following Pigulski et al. (2016). The other stars' data had already been reduced by previous users, though not published, and thus were used without further modification. As an example, we present the MOST lightcurve of the WN7o star WR120 in Figure 3.

Zoom In Zoom Out Reset image size

All the lightcurves analyzed in this work are made available as online material. In order to characterize them, we first carried out period searches to identify the dominant frequencies present in the photometry of each star. For this, we used the discrete Fourier transform method implemented in the PERIOD04 software (Lenz & Breger 2005) to obtain periodograms for each data set of all stars in our sample. These are also available as online material.

As we are interested in the stochastic variability of these WR stars, we first had to remove any a priori known periodicities from our observations before carrying out our analysis, a process known as prewhitening. This was done by evaluating the amplitude signal-to-noise ratio (S/N) of each peak corresponding to a known periodicity in the stars' periodograms. The S/N is calculated as the ratio of the peak amplitude over the average amplitude of residuals inside a symmetric 1 day⁻¹ window around the frequency of the peak, following the methodology given in the Methods section of Bowman et al. (2019a). If the S/N value of a peak is above four, the peak is deemed significant, as suggested by Breger et al. (1993), and the corresponding signal is removed by fitting the sinusoidal model,

$$\begin{eqnarray}&&{\rm{\Delta }}m=A\cos (2\pi \nu (t-{t}_{0})+\phi ),\end{eqnarray} \tag{ 1 }$$

on the lightcurve of time, t, normalized to a zero-point, t₀, with a phase-shift, ϕ, while using the peak's frequency ν and amplitude A.

The BRITE sample contains two stars, WR6 and WR22, with known periodic signals which are expected to exhibit distinct individual peaks in their periodograms. WR6 presents variations in photometry, polarimetry and spectroscopy with a period of 3.76 days that could be explained by the presence of CIRs in its wind (e.g., Morel et al. 1997), although, e.g., Schmutz & Koenigsberger (2019) prefer a binary scenario. In view of the epoch dependency of the variability and the length of the lightcurve, a tailored iterative prewhitening procedure was conducted for WR6 by first carrying out a time-dependent frequency (TF) analysis (also known as Gabor transform) to localise the time intervals over which the various harmonics of the 3.76 day signal are present. We then removed the harmonics with S/N over four by fitting the sinusoidal function previously described in the corresponding time interval and subtracting the fit from the portion of the lightcurve containing the signals. The amplitude and frequencies of the removed harmonics are shown in Table 4. Also, WR22 is a well-known eclipsing binary with a period of 80 days, whose eclipse has recently been modeled by Lenoir-Craig et al. (2021). We used their best-fit Model 1 using the O9V spectral class for the companion to remove two eclipses from the 2017 data and one from the 2018 observation run.

In the TESS sample, eight stars are known to show periodic signals that could be detected within TESS's ∼27 days sectors: WR1, 6, 22, 46, 66, 100, 103, and 134. WR1 and WR134 are both highly likely to harbor CIRs in their winds, with the former exhibiting variations in both spectroscopy and photometry with a period of 16.9 days (Chené & St-Louis 2010), while the spectroscopic variability of the latter shows a clear periodicity of 2.3 days (McCandliss et al. 1994; Morel et al. 1999). The causes of the observed short-term periodic variabilities of WR46 (P ∼ 8 hr) and WR66 (P ∼ 4 hr) are still uncertain, with plausible scenarios ranging from nonradial pulsations to close binaries with compact objects (Antokhin et al. 1995; Hénault-Brunet et al. 2011; Gosset et al. 2011; Skinner et al. 2021). Spectroscopic observations of WR100 were conducted by Chené & St-Louis (2011), and CIR-type large-scale line-profile variability has been observed in its winds, albeit without the determination of a clear period. Chené & St-Louis (2011) also investigated spectroscopic observations of WR103, and small-scale line-profile variability was detected, again with no clear period detected. Except for WR103, significant peaks were detected near the expected frequencies (and sometimes harmonics) for all of these stars, with the frequencies and amplitudes of the removed signals shown in Table 4. Also, we completely removed the eclipse data from the TESS sector 36 observations of WR22, as the shape of the TESS-observed eclipse is sufficiently different from the ones observed by BRITE to require it to be numerically modeled separately to ensure a correct detrending, which goes beyond the scope of the present study.

Three stars in the MOST sample are also known to exhibit periodic signals: WR103, 115, and 120, with the latter two exhibiting large-scale CIR-type variability in spectroscopic observations (Chené & St-Louis 2011), with no clear detected periodicity. No significant peak were detected in the periodograms of any of the corresponding data sets.

Once the known significant periodicities were removed, the same procedure was iteratively repeated for all other peaks with S/N bigger than four in all of our sample stars' periodograms, again following the details given in the Methods section of Bowman et al. (2019a). Many more significant signals were removed from multiple TESS lightcurves, with their frequencies and amplitudes shown in Table 4.

A thorough analysis of all the WR periodograms in our sample revealed no other outstanding peak than those described in Table 4. The periodograms all show a forest of a large number of peaks with an increasing amplitude toward lower frequencies, starting from ∼0.5 day⁻¹ down to ∼0.1 day⁻¹, sometimes followed by a decrease in amplitude below 0.1 day⁻¹, similarly to what had been observed by Ramiaramanantsoa et al. (2019) for the WN8 star WR40. This decrease in amplitude is presumably intrinsic to the stars, being unrelated to any details of the detrending process. An example of an amplitude spectrum we obtained for the WN6ha star WR24 is shown in Figure 4. All amplitude spectra are available as online material.

Zoom In Zoom Out Reset image size

In order to characterize the stochastic, low-frequency variability for our sample of WR stars, we fitted a combination of red and white noise to all our amplitude spectra as described, for example, in Bowman et al. (2019a) or Bowman et al. (2019b), by applying a Markov Chain Monte Carlo algorithm using the Python code Emcee (Foreman-Mackey et al. 2013) to fit the amplitude spectrum of each lightcurve with a semi-Lorentzian distribution:

$$\begin{eqnarray}&&\alpha (\nu )=\displaystyle \frac{{\alpha }_{0}}{1+{\left(\tfrac{\nu }{{\nu }_{\mathrm{char}}}\right)}^{\gamma }}+{C}_{W},\end{eqnarray} \tag{ 2 }$$

where α₀ represents the amplitude of the frequency-dependent component in the limit of zero frequency; ν_char is the characteristic frequency of the stochastic variability present in the lightcurve, which is the inverse of the characteristic timescale following ν_char = (2π τ)⁻¹; and γ is the logarithmic amplitude gradient. The first term on the right-hand side of this equation corresponds to the frequency-dependent red noise, while the second term, C_W , is the frequency-independent white noise. It is important to note that the values of fitted α₀ are significantly higher than the error bars of the photometric observations, and therefore are actually genuine intrinsic signals from the star and not noise in the traditional sense (see Bowman et al. 2019a). The white-noise components, however, are of much smaller amplitude, lower than or of the same level as the error bars, and are more likely to correspond mainly to instrumental or photon noise. Similarly to Bowman et al. (2019a), we then tested the goodness of the fit against a white-noise-only model using statistical F-tests. In all cases, we found the red-noise component to be significant. An example of such a fit to an amplitude spectrum is shown in Figure 4 for WR24 superposed on the observations plotted in black. The red- and white-noise components are plotted as red and blue dashed curves, respectively, while the resulting fit is shown as a solid green curve. All other fits can be found superposed on the amplitude spectra in the online material, and all the fitted parameters are listed in Tables 1, 2 and 3. We note that the values of our fitted parameters are in agreement with the range of red-noise parameters for WR stars shown in Figure 4 of the recent work done by Nazé et al. (2021) on the presence of red noise in the photometric data of evolved massive stars.

In order to check for trends between the values of the fitted model parameters and multiple stellar characteristics, we first had to ensure the comparability of the data. The data we are using in this study have been acquired by three different telescopes observing with various cadences in neighboring and/or overlapping wavelength ranges in the visible and very-close (700–1000 nm) infrared parts of the spectrum over different observation periods. Fortunately, it is well established that the source of visible continuum light for WR stars is wavelength-independent electron scattering in the winds, and therefore we do not expect that the different observation bands will significantly impact the observed fluxes, periodograms and fitted model parameters. Also, the duration of the observational data determines the lowest detectable frequency in the amplitude spectra, which for data covering multiple weeks will always be below the frequency range over which the red noise is observed (∼0.1–0.5 day⁻¹). Therefore, this should also not significantly affect the fitted parameters.

However, the different sampling rates of the telescopes are expected to influence the values of the fitted model parameters, and a vital step here is to quantify by how much. To do this, we used all TESS lightcurves with a 2 minute cadence in our sample and created as many lightcurves with a uniform sampling of 10, 30, and 100 minutes as possible from them. We had to exclude the 2 minute lightcurves from WR3, 4, 33, and 60 from the 100 minute lightcurve creation process, as some of the lightcurves created from them showed no red noise in their periodograms. We were thus able to create 170, 510, and 1500 lightcurves with, respectively, 10, 30, and 100 minute cadences from our 34 original 2 minute lightcurves. We then fitted the red and white noise semi-Lorentizan model from Equation (2) on their periodograms, extracted the optimal model parameters, and normalized α₀, ν_char, and γ by the original best-fit parameter value of the 2 minute lightcurve. We show in Figure 5 the averages and standard deviations of all parameters for the different cadences.

Zoom In Zoom Out Reset image size

We observe from Figure 5 that the fitted value of α₀ seems to follow a negative trend with increasing standard deviation as the cadence goes down. The 10 minute lightcurves have ${\alpha }_{0}\sim (0.96\pm 0.06){\alpha }_{2\min }$, the 30 minute lightcurves have ${\alpha }_{0}\sim (0.90\pm 0.12){\alpha }_{2\min },$ and the 100 minute lightcurves have ${\alpha }_{0}\sim (0.83\pm 0.22){\alpha }_{2\min }$. We remark that the means of the 10, 30, and 100 minute cadences always stay within one standard deviation of the 2 minute value. From this, we can estimate an intrinsic uncertainty on the value of α₀ for cadences lower than 2 minutes, which is ±6% for 10 minute cadences, ±13% for 30 minute cadences, and ±27% for 100 minute cadences. However, as shown in Figures 6 and 7, these uncertainties on α₀ for lightcurves with 10, 30, and 100 minute cadences are negligible when compared to the uncertainty brought forth by the varying values of fitted parameters over different epochs of observation, as the structures in the winds of the WR stars are changing over time, affecting the level of electron scattering.

Zoom In Zoom Out Reset image size

ν_char appears to be the parameter least influenced by a changing cadence in Figure 5, with a relatively constant mean and increasing uncertainty as cadence goes down. It is important to note that the error bars of the average 100 minutes value of ν_char, γ, and C_W are overestimated due to their non-normal distributions. The γ and C_W parameters exhibit a similar behavior to one another, both showing an increase in their values and spread as the cadence gets lower. This was expected for C_W , since white noise is known to be dependent on the satellite cadence of observations and on the photometric precision (Bowman et al. 2020). As with α₀, the values of the intrinsic errors on the parameters due to the cadence for the 10, 30, and 100 minute lightcurves are negligible when compared with the typical variations of these parameters between different epochs, as shown in Figures 6 and 7.

Zoom In Zoom Out Reset image size

We then looked for correlations between the parameters describing the red-noise component of the photometric changes observed in our stars with various stellar and wind parameters. In order to do so, we first separated our sample into two separate groups, one including stars that showed signs of the presence of hydrogen in their spectra, and the other including stars that did not. In Figures 6 and 7, we present log-log plots of α₀, ν_char, and γ, all provided in Tables 1–3, as a function of the wind terminal velocity, mass-loss rate, stellar temperature, absolute magnitude in the Smith v band, ${{\rm{M}}}_{v}^{\mathrm{WR}}$, stellar radius, and bolometric luminosity. The stellar radius, R_*, in this context refers to the location where the Rosseland continuum optical depth is equal to τ_Ross,c = 20, and the stellar temperature, T_*, is the fictitious quantity calculated using the Stefan–Boltzmann relation at that position, as was done in Hamann et al. (2019) and Sander et al. (2019). TESS data are in blue, BRITE data in black and MOST data are in red, with WC stars indicated by circles, WN stars by inverted triangles and the WO2 star WR142 by a square. The wind and stellar parameters were taken from Hamann et al. (2019) for WN stars and Sander et al. (2019) for WC and WO stars, and the absolute magnitudes were taken from Rate & Crowther (2020). For WR79b, the stellar parameters were taken from Bohannan & Crowther (1999), and for WR138a from Gvaramadze et al. (2009), as these stars did not figure in the previously mentioned references. For the three stars in our sample for which no individual values were available, namely WR8, 87, and 159, we used the average value for stars of the same spectral type from Hamann et al. (2019). Note that recently there have been great strides made in the understanding of how the optically thick winds of WR stars are launched and driven. Grassitelli et al. (2018) presented one-dimensional hydrodynamic stellar structure models for He stars that demonstrate that outflows can be driven to supersonic velocities by the opacity associated to the recombination of Fe atoms, the so-called Fe-peak opacities (at T ∼ 200kK) and/or the He ii opacity bump (at T ∼ 50 kK). They find that above a given $\dot{M}$, the winds can be driven entirely by Fe. For winds that cannot, the additional He ii opacities result in an extended low-density envelope accompanied by a density inversion. The stars that launch these Fe-driven winds all have a very compact core, close to 1 R_⊙. Until very recently, this result was in contrast with results of stellar atmosphere codes that estimated radii higher by a factor of 2–3. However, Sander et al. (2020) have now produced improved atmosphere models for classical WR stars in which the hydrodynamical equation of motion is solved throughout the atmosphere. They find, in accord with Grassitelli et al. (2018), that the winds are driven by the peak in the Fe opacities. Furthermore, they also show that the complex ionization structure of the wind that helps to drive the wind in the outer layers results in complex velocity laws that cannot be well described by a simple beta law. Similar conclusions were reached by Poniatowski et al. (2021) using a combination of static OPAL opacities usually used in static stellar structure models with a simple parameterization equation to take into account the additional line opacity from a supersonic flow. These velocity laws are found to include bumps, plateaus, and even deceleration regions. This is an important result because the stellar parameters, and in particular the stellar temperature, that have been derived in the papers we have used are based on using a beta law to relate quantities in the observable part of the wind to those at the stellar surface. At the present time, no updated values of these stellar parameters are available and it is unclear what the effect of assuming a beta law is on the published parameters ¹¹ we have used to search for correlations. Obviously, once new values become available, the study will have to be redone.

We fitted a straight line to all these data sets using a least-squares algorithm and overplotted the result as a dashed black line. The value of the slope and its corresponding error is given in the top-right corner of each panel along with the Spearman's rank R_s coefficient. If more than one data set existed for a given star, we linked the fitted values of the parameters describing the red-noise component with a vertical line to show the range covered by the values we obtain for the various observations.

There are seven graphs in Figures 6 and 7 that have slopes above the 3σ level and a two-sided Spearman's rank coefficient that is significant at p < 0.001 for the number of points in our samples. However, only the negative trend between $\mathrm{log}{\alpha }_{0}$ (we refer to log₁₀ as $\mathrm{log}$ throughout the text) and logv_∞ also has a significant p-value at the p < 0.05 level associated with an F-test between the χ² of the fitted linear trend and the χ² of a best-fit flat model for both cWR and WNh stars. Thus, although we can determine the presence of correlations for the six other graphs, we cannot ascertain their statistical significance with our current sample.

Correlations satisfying the Spearman's rank p < 0.001 criterion for the correlation coefficient, r_s , are observed between $\mathrm{log}{\alpha }_{0}$ and both $\mathrm{log}\dot{{\text{}}{\rm{M}}}$ and $\mathrm{log}\,{{\rm{L}}}_{\mathrm{bol}}$ in Figure 6, indicating that hotter (since v_∞ and T_* are related), more-luminous cWRs with denser winds appear to exhibit lower levels of stochastic variability than their cooler counterparts. We note here that the Spearman's rank of $\mathrm{log}{\alpha }_{0}$ versus $\mathrm{log}\dot{{\text{}}{\rm{M}}}$ is just slightly under the threshold. We also note that the negative correlation between $\mathrm{log}{\alpha }_{0}$ and $\mathrm{log}\,{\text{}}{{\rm{L}}}_{\mathrm{bol}}$ is contrary to the results of Bowman et al. (2020) for OB stars, who found that higher-luminosity stars tend to be characterized by larger values of α₀.

For WNh stars, correlations are observed between $\mathrm{log}{\alpha }_{0}$ and $\mathrm{log}{\text{}}{{\rm{T}}}_{* }$, and also between $\mathrm{log}{\nu }_{\mathrm{char}}$ and both $\mathrm{log}{\text{}}{{\rm{v}}}_{\infty }$ and $\mathrm{log}{\text{}}{{\rm{T}}}_{* }$. An important remark to make here is that these last two trends drop below of the Spearman's rank coefficient if the TESS data from WR3 is removed from the sample, as its two data sets observed over different periods both have a fitted γ > 9, by far the highest values of our sample, with the next highest value being γ ∼ 4.7, and the overall (cWR + WNh) mean and standard deviation being 1.7 ± 0.8. Since a close inspection of the data and procedure has not revealed anything anomalous, and since the periodograms of both independently observed epochs have been fitted with such a high γ, we consider it real and not an erroneous observation. Note also that T_* and v_∞ are known to be correlated: hotter stars tend to have faster winds, although this has not been formally quantified (e.g., Hamann et al. 2019). As previously mentioned, we can also observe in Figures 6 and 7 that for stars that have been observed over different epochs, the fitted coefficients can be significantly different from one observation to the next, even when they are observed with the same satellite. We consider that this is real scatter and most likely represents a change in the conditions in the wind that must be taken into account. We therefore estimated a standard deviation characterizing these changes as follows:

$$\begin{eqnarray}&&\overline{\sigma }=\displaystyle \frac{\sum {n}_{i}{\sigma }_{i}}{N},\end{eqnarray} \tag{ 3 }$$

where n_i is the number of different observations of the ith star, ${\sigma }_{i}=\sqrt{\sum \tfrac{{(\overline{x}-{x}_{i})}^{2}}{{n}_{i}-1}}$ corresponds to the standard deviation of the chosen parameter calculated over the n_i observations of that same star, while N corresponds to the total number of observations of stars that have been observed multiple times. The resulting values are plotted in the bottom-left corner of each panel in Figures 6 and 7. Using these values as error bars for our F-tests, we find one significant correlation at p < 0.05 between α₀ and v_∞ for all WR stars, with corresponding p-values of 0.02 for cWR and 7 × 10⁻⁴ for WNh. We believe that this trend should be investigated further, as should be the other observed correlations once more data sets are acquired and the range of intrinsic change of the various parameters is better established. We also note that a similar trend between spectroscopic variability and T_* was found in spectroscopic data sets of WR stars by Chené et al. (2020).

In order to obtain a different perspective on the timescales of the variability present in our photometric observations, we carried out wavelet analyses on all data sets for the stars in our sample following the method of Torrence & Compo (1998). An example is shown in Figure 8 for the TESS data of WN6ha star WR24. The top panel shows a reconstruction of the lightcurve of this star using 60 minute bins, and the bottom panel shows the wavelet power-spectrum of this data set. We find that the stronger signals have periods in the 1–5 days range. This is consistent with the results of our amplitude spectrum analysis presented in the previous section, for which the average value of the characteristic frequencies, ν_char, for the red-noise component of the variability of cWR stars is 1.2 ± 0.9, while for WNh stars this range is 1.7 ± 1.5 (see middle panels of Figures 6 and 7).

Zoom In Zoom Out Reset image size

As the LDI is normally assumed to be intrinsic to the driving of hot-star winds, could this process alone explain the observed photometric variability we found in WR stars? And, if not, can we determine if the LDI contributes significantly to the variability?

Sundqvist et al. (2018) have shown that for a wind with a velocity structure described by a beta velocity law, the LDI produces clump sizes of the order of the Sobolev length, ${L}_{\mathrm{Sob}}={v}_{\mathrm{th}}/({dv}/{dr})\propto {v}_{\mathrm{th}}/({v}_{\infty }/{R}_{* })\propto {R}_{* }\sqrt{{T}_{\mathrm{wind}}}/{v}_{\infty }$. Note that although here the thermal speed is related to the wind temperature and not the stellar temperature, stellar and wind temperatures are correlated. Since the thermal speed is expected to be about 1% of the terminal speed, this process is expected to lead only to very small clumps, according to the above equation. It is far from clear if such small structures can lead to the level of photometric variability we see in WR winds.

On the other hand, if WR winds are self-similar, a more appropriate parameter to consider would then be L_Sob/R_*, which is proportional to the square root of the wind temperature divided by the terminal velocity. As faster winds are known to be hotter (Sander et al. 2019), it is possible, if the correlation between v_∞ and T_* is steeper than the square root, that the Sobolev length inversely correlates with the temperature, i.e., hotter winds have smaller Sobolev length and therefore potentially smaller wind structures. This in turn can lead to a lower level of variability as larger clumps are expected to scatter more flux, and while a high number of smaller clumps can possibly scatter the same amount of flux, they would need to occupy a larger portion of the surrounding area of the WR disk to do so, a configuration we think is less likely to happen than single, larger clumps doing the scattering.

In order to check that increasing clump size directly leads to increased variability, we used the same simple clumping model as presented in Ramiaramanantsoa et al. (2019). These authors where able to produce similar lightcurves to that observed with BRITE for the WN8h star WR40 by populating its wind with a distribution of clumps and calculating the resulting continuum light that is Thomson scattered off these clumps as they travel through the wind. The f parameter in this model sets the percentage of incoming stellar flux that is Thomson scattered by the biggest clump. Smaller clumps scatter less flux, with their size following a negative-index power law with exponent −2/3. Using increasing values of f to simulate the increasing clump size, multiple synthetic lightcurves were produced and their amplitude spectra were extracted. Following the same recipe as with our sample stars, our artificial lightcurves were then fitted with the semi-Lorentzian distribution previously shown. Figure 9 shows the fitted amplitude of variability plotted against the various f values used.

Zoom In Zoom Out Reset image size

A clear positive trend appears, showing that, according to this simple model, stars with bigger clumps in their winds are expected to show an increased amplitude of variability in their amplitude spectrum, similar to the observed trend.

In Figure 10 we present the values found in our simulations for the steepness parameter of the red component of the stochastic variability, γ, as a function of various model parameters: f; the lifetime of clumps, T; the mass-loss rate, $\dot{M};$ the stellar radius, R_*; and the exponent of the size distribution, $\exp$. Although no obvious correlation is found, we note that all the lightcurves produced by the model have an associated γ value in the range 1.5 < γ < 3.0, with the exception of simulations with a clump lifetime parameter T below 2 hr, for which the range becomes 3.5 < γ < 5.0. As can be seem from Figures 6 and 7, this partly overlaps with most of the γ values determined from our data sets, which are concentrated in the range 0.8 < γ < 3.0, with a smaller number with γ > 3.0.

Zoom In Zoom Out Reset image size

Krtička & Feldmeier (2018, 2021) have recently suggested that variable wind blanketing that results in extra heating of the stellar photosphere can explain the stochastic variability in optically thin winds of O stars. Indeed, a variable mass-loss rate caused by LDI can in turn produce detectable light variability at the photosphere of the order of tens of millimagnitudes. Although for WR stars this effect is unlikely to contribute directly to the photometric variability since the stellar surface is hidden by the optically thick stellar outflow, it might propagate in the wind and contribute somehow to the wind variability in the observable parts of the wind.

Therefore, there are several arguments in favor of LDI playing an important role in generating the observed stochastic photometric variability. LDI could even lead to an anticorrelation between the level of the photometric and spectroscopic stochastic variability, as observed in the winds of WR stars, with hotter winds having smaller clumps, leading to a lower amplitude of variability. However, this does not exclude that an additional process occurring at the stellar surface could seed the instabilities and affect the clump characteristics in the winds. In the next two sections, we discuss two such possible processes.

Following the original work of Cantiello et al. (2009) on the effects of subsurface convection zones caused by an FeCZ on the optically thin winds of OB stars, Grassitelli et al. (2016) studied the effects of this convective zone on the radiative envelopes of WR stars in the 2–17 M_⊙ mass range. Although only small-amplitude (∼1 km s⁻¹) total velocity perturbations were found for lower-mass stars, they found that for masses larger than ∼10 M_⊙, surface convective velocities of the order of 10 km s⁻¹ were expected. More precisely, they predicted a linear increase in convective velocities with stellar mass from 2 to 26 km s⁻¹, approaching the local sound speed for the higher masses. Although these turbulent surface velocity fields are predicted to be somewhat attenuated in the radiative zones above the FeCZ, the predicted surface total velocities, v_surf, are of the same order of magnitude as the convective velocities and are predicted to follow the same increase with mass. These predictions are supported by observations of stochastic spectroscopic variations of the winds of WN stars, quantified by the rms of the line flux relative to the continuum, that are found to show an increase in amplitude with increasing mass (see their Figure 10).

The range of dominant frequencies of the variability (ν_char) can help shed more light on the physical mechanism behind the observed variability. FeCZ are expected to have a typical turnover frequency of 6–60 μHz (0.5–5 day⁻¹) (Cantiello et al. 2009; Cantiello & Braithwaite 2011, 2019), which can be transmitted to the above radiative layer since it is below the Brunt–Väisälä frequency of 0.1–1 mHz (Cantiello et al. 2009). As shown in Figures 6 and 7, almost all the stars in our sample have a characteristic frequency in this expected range. However, the typical timescale associated with convection determined for helium stars, 10²–10³ s (Grassitelli et al. 2016) is some two orders of magnitude shorter, leading to much higher characteristic frequencies of $\mathrm{log}{\nu }_{{conv}}=2\mbox{--}3\,{\mathrm{day}}^{-1}$. This does not agree with our observed values.

The temperature of the outer hydrostatic layers of OB and WR stars is expected to be an important parameter in the surface variability produced by the FeCZ (Cantiello et al. 2009; Grassitelli et al. 2015) because it has a direct influence on the location of the FeCZ caused by the partial ionization of Fe-peak elements at $\mathrm{log}{\text{}}{\rm{T}}\approx 5.3$. For higher temperatures, the FeCZ will be located closer to the surface, and therefore the observed photometric variability will be reduced because the convective mass flux perturbing the base of the wind will be smaller. This effect could play a role in explaining the origin of the possible trend we find in our sample between the amplitude of photometric variability and the stellar temperature, T_eff, if confirmed (see Figures 6 and 7). This supports the idea that a subsurface convection zone can directly influence the formation of clumps in the winds of WR stars, particularly for stars above a mass of ∼10 M_⊙.

Another possible source of variability in hot, massive stars that has been put forward in the recent literature is a large number of IGWs, excited by turbulent convection at the interface of the convective core and the radiative zone lying above it (Aerts & Rogers 2015; Bowman et al. 2019a, 2020). IGWs have been claimed to be one possible source of waves that can produce the velocity field required to generate the ubiquitous macroturbulent broadening that characterizes the absorption lines of many OB stars, with the presence of a subsurface convection zone being another (e.g., Aerts & Rogers 2015). In support of this interpretation, a strong correlation was found between the surface macroturbulent speed and the amplitude of stochastic variability in OB stars; see Bowman et al. (2020), who also showed that IGWs become increasingly important in stars with larger spectroscopic masses and luminosities. Note, however, that in order to avoid numerical instabilities, the thermal diffusivity must be increased in those models, and as a consequence the heat flux must be artificially increased to compensate for the resulting damping (e.g., Rogers et al. 2013; Aerts & Rogers 2015). Other authors also claim that the observed photometric variability of OB stars is incompatible with the predicted signature (e.g., Shiode et al. 2013; Lecoanet et al. 2019). As WR stars are among the most luminous stars, IGWs could potentially play a role in their internal structure, as well. However, contrary to massive MS stars, we do not see the surface of WR stars directly, and therefore we can only see these perturbations if they propagate outward in the wind. For the optically thin wind of the O9.5 supergiant star HD188209, Aerts et al. (2017) have shown that perturbations formed at the stellar surface do indeed seem to propagate in the wind. Similarly, this was also demonstrated for the B1 Iab supergiant ρ Leo (Aerts et al. 2018). It remains to be shown what happens in the optically thick winds of WR stars.

The steepness of the amplitude spectrum of the stochastic photometric variability (characterized by the γ exponent) can also help shed some light on the cause of the stochastic variability. Recent three-dimensional simulations by Edelmann et al. (2019) predict that IGWs are compatible with a steepness of amplitude exponent 0.8 < γ < 3.0, with the value also depending on the rotation rate and evolutionary stage of the star (Ratnasingam et al. 2020). As can be seen from Figures 6 and 7, we find that the vast majority of WR stars in our sample have a γ exponent in that range.