© The Authors 2025

1. Introduction

Massive stars are progenitors of neutron stars and black holes, and play a central role in galactic chemical evolution and star formation (Kippenhahn & Weigert 1990; Bromm & Larson 2004). Constraining their evolution is key to understanding both compact object populations and galaxy evolution. However, several critical processes remain poorly understood. Rotation, in particular, influences stellar evolution through internal mixing and angular momentum transport (Maeder & Meynet 2000; Langer 2012; Aerts et al. 2019), but its measurement via spectroscopy is complicated by an additional source of broadening known as macroturbulence.

Macroturbulence is often attributed to large-scale velocity fields in stellar atmospheres, though its physical origin remains poorly understood and unconstrained (Aerts et al. 2009; Grassitelli et al. 2015; Simón-Díaz et al. 2017). Three main hypotheses have been proposed to explain the high macroturbulent velocities observed in massive stars. Aerts et al. (2009) suggested that the collective effect of numerous low-amplitude gravity-mode pulsations can give rise to the additional broadening of the spectral line profile, which is usually attributed to macroturbulence. Grassitelli et al. (2015) identified a correlation between macroturbulence and turbulent pressure derived from 1D stellar structure models, proposing that pressure fluctuations in sub-surface convection zones can drive small-amplitude pulsations, which in turn result in the observed velocity fields. More recently, Debnath et al. (2024), using 2D hydrodynamic simulations, demonstrated that gas pockets near the iron-opacity bump can be radiatively accelerated to velocities comparable to those needed for macroturbulence. Previous spectroscopic studies have suggested that the time-averaged effects of pulsations generates the velocity field responsible for macroturbulence (v_macro) in massive stars (Aerts et al. 2009, 2014; Aerts & Rogers 2015). Moreover, massive stars typically have large values of v_macro and are largely located within theoretical pulsation instability regions in the Hertzsprung–Russell diagram (Simón-Díaz et al. 2017; Godart et al. 2017; Burssens et al. 2020). However, the time-dependent relationship between macroturbulence and pulsations, as well as the impact on inferred atmospheric parameters, remains unclear.

To date, only a handful of studies have investigated line profile variability (LPV) in massive stars and its connection to multi-period pulsations and macroturbulence (see Bowman 2020 for a review). One example is the nearest O-type star, ζ Oph (Kambe et al. 1997; Walker et al. 2005; Howarth et al. 2014). Time-series high-resolution spectroscopy from various campaigns has allowed a dominant period of 3.33 h (i.e. ν = 7.18 d⁻¹) to be consistently observed and attributed to a non-radial pulsation mode with spherical harmonic indices ℓ = |m| = 4 (Kambe et al. 1990; Reid et al. 1993; Kambe et al. 1997; Balona & Kambe 1999). Other reported pulsation periods include 2.01 h (ν = 11.89 d⁻¹, ℓ = 8; Balona & Kambe 1999; Kambe et al. 1997), 5.34 h (Kambe et al. 1997), 4.468 h (ν = 5.37 d⁻¹), and 2.43 h (ν = 9.86 d⁻¹; Reid et al. 1993). However, none of these investigations explored how such pulsations are related to the time-dependent nature of macroturbulence.

Space-based photometry from the Microvariability and Oscillations of Stars (MOST; Walker et al. 2003) satellite further confirmed the multi-period pulsations of ζ Oph, detecting a persistent signal with a period of 4.63 h (i.e. ν = 5.18 d⁻¹) across several years (Walker et al. 2005). Using SMEI (Solar Mass Ejection Imager; Eyles et al. 2003), WIRE (Wide-field Infrared Explorer; Schember et al. 1996), and MOST space photometry, Howarth et al. (2014) corroborated the detection of the 4.63 h and 3.33 h periods and revealed additional pulsations with periods of 8.11 h (ν = 2.96 d⁻¹) and 8.99 h (ν = 2.67 d⁻¹). Therefore, it is well established that ζ Oph is a multi-period pulsator, and that its variability is not strictly periodic (Howarth et al. 2014). Moreover, the connection between its photometric and spectroscopic pulsation periods warrants further investigation.

In this paper we analyse the pulsations and various diagnostics of LPV, including macroturbulence, for ζ Oph using new TESS (Transiting Exoplanet Survey Sattelite; Ricker et al. 2015) photometry and high-resolution spectroscopy from the European Southern Observatory (ESO) ESPRESSO spectrograph (Pepe et al. 2021). ζ Oph is a rapidly rotating, pulsating O9.2 IVnn star (Sota et al. 2014) with an uncertain projected rotational velocity of 320 < v sin i < 400 km s⁻¹ (Conti & Ebbets 1977; Herrero et al. 1992; Jankov et al. 2000; Simón-Díaz & Herrero 2014), which corresponds to a rotation period of between 0.7 and 1.2 d. This uncertainty arises primarily because sparsely sampled spectra cannot disentangle the time-dependent contributions of rotation, pulsations, and macroturbulence in the line profiles. With unprecedented temporal coverage, we empirically measured the amplitudes and timescales in the line-specific variability in radial velocity (RV), rotation rates, and macroturbulence induced by pulsations. We also investigated the impact of variable macroturbulence on the inferred stellar parameters, and potential biases in reported values of macroturbulence from the longer integration times needed for massive stars fainter than ζ Oph.

2. Data description and reduction

Our spectroscopic data were obtained with the ESPRESSO spectrograph (Mégevand et al. 2014; Pepe et al. 2021) at the combined Coudé facility of the Very Large Telescope (VLT) of ESO in Chile, under programme 113.26B9 (PI: Bowman). All observations were conducted in the Ultra High-Resolution 1-UT mode, providing a median resolving power of 190 000 and a wavelength coverage of 380–788 nm. This configuration ensures a RV precision of ∼5 m s⁻¹. A total of 1706 spectra were assembled, each with a 30-s exposure time over four consecutive nights of 7 − 10 June 2024. The ESPRESSO spectra were reduced using the ESOReflex pipeline¹ (Freudling et al. 2013). The reduced spectra were manually normalised using spline fitting to the continuum level. The mean signal-to-noise ratio (S/N) across all spectra at the wavelength of 550 nm is 393, with the 16th and 84th quantiles for S/N being 295 and 495, respectively.

The ESPRESSO spectra of ζ Oph are rich in metal lines, but given its rapid rotation rate (Herrero et al. 1992; Howarth & Smith 2001; Simón-Díaz & Herrero 2014) most metal lines exhibit significant blending. However, the He I+ IIλ4026, He Iλ4143, He IIλ4200, He Iλ4471, and He Iλ4922 are sufficiently isolated to perform a LPV study unaffected by blending. The Balmer lines are strongly affected by Stark broadening and are less sensitive to non-radial pulsations than helium lines. For these reasons, it is typical to focus on isolated helium lines in LPV studies of rapidly rotating massive stars (e.g. Balona & Kambe 1999). In this paper we show results for the He I+ IIλ4026, He IIλ4200, and He Iλ4922 lines and investigate how rotation, pulsations, and macroturbulence vary among different ionisation stages of helium. We leave the analysis of blended metal lines for future work.

3. Line profile variability analysis

Our ESPRESSO time-series of ζ Oph reveals significant LPV caused by multi-period pulsations. In Fig. 1 we show three spectra selected from 9 June 2024 to demonstrate the variability in three helium lines, as well as the mean spectrum obtained by averaging all 1706 spectra overplotted in dark blue for each line. Figure 1 also shows residual dynamic spectra, constructed by subtracting the mean spectrum from each individual spectrum. The residual spectra are time-stacked along the y-axis and reveal LPV (i.e. ‘bumps’ moving across the line profile) caused by non-radial pulsations (see e.g. Aerts et al. 2009). In the following, we demonstrate how pulsation-induced LPV translates into variability in other spectroscopic diagnostics.

Fig. 1. Left panels: three spectra from 9 June 2024 showing LPV for the He I+IIλ4026, He IIλ4200, and He Iλ4922 lines (from top to bottom), with the stacked mean spectrum overplotted in dark blue. Right panels: residual dynamic spectra for the corresponding left panels.

3.1. Radial velocity

The RV stability of individual spectral lines was examined using the cross-correlation technique (Zucker 2003; Shenar et al. 2020). Figure 2 shows the observed RV variability for the same spectral lines shown in Fig. 1, which all exhibit peak-to-peak RV variations exceeding 20 km s⁻¹. We investigated the dominant variability periods in the RV time series using an iterative pre-whitening technique (see Bowman & Michielsen 2021). This method involves fitting a sinusoidal model of the form

$$\begin{array}{c}\hfill f(t)={\displaystyle \sum _{i=1}^n}{A}_{i}sin({\omega }_i(t-t_0)+{\varphi }_i),\end{array}$$(1)

Fig. 2. RV time series for three helium lines in ζ Oph extracted using cross-correlation. The peak-to-peak RV variations exceed 20 km s−1 in a multi-period manner because of pulsations rather than binarity.

where A_i is the amplitude, ω_i is the frequency identified from the Lomb-Scargle periodogram analysis of the time series (see Scargle 1982), t is the time, and ϕ_i is the phase relative to the reference time t₀, chosen as MJD = 60468.983772, which is the time stamp of the first spectrum. The pre-whitening process was continued until the S/N of the next highest peak in the frequency spectrum dropped below a false alarm probability of 0.001%. At each step, we carefully examined the corresponding spectral window of the ESPRESSO time series (see Fig. 8) to mitigate aliasing. After all significant periods had been identified and extracted, they were included in the model in Eq. (1) and optimised using least-squares fitting. This method is typical for analysing light curves of pulsating massive stars (see the review by Bowman 2020).

Fig. 8. Spectral windows for the ESPRESSO spectroscopic time series (top panel) and the TESS photometric data (bottom panel) analysed in this work.

Iterative pre-whitening of the RV times series of ζ Oph revealed a range of variability timescales and amplitudes, some of which were reported in previous studies. We provide the significant periods and their corresponding uncertainties from the RV (and other) time series analysis in this work in Table B.1. The 3.33 h period was observed for all three helium lines, which are shown in Fig. 2. In the He IIλ4200 and He Iλ4922 lines, the 3.33 h period was the dominant source of variability with an amplitude of 3 km s⁻¹, while for the He I+ IIλ4026 line, this was the second largest source of variability with an amplitude of 2 km s⁻¹. This period was previously shown to be due to a non-radial pulsation mode with spherical harmonic geometry ℓ = |m| = 4 (Kambe et al. 1990; Reid et al. 1993; Kambe et al. 1997). We also find a significant period of 2.04 h in both the He I + IIλ4026 and He Iλ4922 lines, while the He IIλ4200 line exhibits a similar period of 2.00 h. Additionally, we detect periods of 2.46 h and 2.50 h in the He IIλ4026 and He Iλ4200 lines, respectively, which are close to the reported 2.43 h pulsation period measured by Reid et al. (1993) and Kambe et al. (1993). Later, Kambe et al. (1997) argued that the 2.43 h period could be an alias of the true pulsation period of 2.01 h. In our ESPRESSO RV time series data, both of these periods are significant. However, in the TESS light curve of ζ Oph, we only find the 2.43 h period to be significant, while the 2.01 h period is below the significance threshold of 0.001% false alarm probability (see Sect. 5). Thus, we conclude that the 2.43 h period is likely the real pulsation mode and the 2.01 h period is an alias, as concluded by Kambe et al. (1997).

It is common for massive stars to be found in binaries, and the majority of them to interact with a companion during their lifetime (Sana et al. 2012; de Mink et al. 2013; Marchant & Bodensteiner 2024). Therefore, we also searched for evidence of binarity in our RV time series. The high surface abundance of nitrogen and helium, along with the rapid surface rotation suggest that ζ Oph is a binary interaction product (Herrero et al. 1992; Villamariz & Herrero 2005). Additionally, ζ Oph is also classified as a runaway star because of its high peculiar spatial velocity of 30 km s⁻¹ (Blaauw 1961), suggesting that it could have been ejected following the supernova explosion of its former primary companion. Based on this premise, Renzo & Götberg (2021) show that an accretor star with an initial mass of 17 M_⊙ in a 100-d binary orbit with a 25-M_⊙ donor can broadly reproduce the observed mass, radius, runaway velocity, surface rotation and surface abundances of helium, nitrogen and carbon of ζ Oph. Therefore, investigating any potential signatures of binarity for ζ Oph remains an important complementary goal of our work.

In most spectroscopic studies of massive stars that focus on binarity detection, pulsations are often treated as a noise term in the RV time series. A peak-to-peak RV amplitude exceeding 20 km s⁻¹ is typically interpreted as evidence of binarity (Sana et al. 2013; Dunstall et al. 2011; Mahy et al. 2022; Bodensteiner et al. 2025). However, we find significant RV variability above this threshold that is clearly driven by pulsations. Moreover, the variability is line-dependent and multi-period, consistent with earlier findings (Kambe et al. 1990, 1997; Walker et al. 2005; Howarth et al. 2014). Our results confirm the quasi-periodic nature of the pulsations: the recovered periods are close to, but not exactly the same as those previously reported. We therefore conclude that the observed RV variability in ζ Oph solely arises from pulsations, as we find no spectroscopic evidence of binarity in the RV time series.

Of course, given the observing cadence and time span of our ESPRESSO time series, we are most sensitive to the presence of a close binary companion. Using the mass of ζ Oph reported in literature (20.2 M_⊙; Repolust et al. 2004) and its observed RV semi-amplitude from our ESPRESSO data (i.e. ∼15 km s⁻¹), we should be able to detect signatures of any binary companion of mass roughly 1 M_⊙ or larger from our RV time series. However, since the total observing window is only 4 d, we can only exclude the possibility of companions with periods of less than about 10 d. Finally, despite the lack of binarity in ζ Oph, our results demonstrate the importance of accounting for LPV caused by pulsations in multiplicity studies of massive stars (see Simón-Díaz et al. 2024).

3.2. Rotation

Fig. 3. Left panel: Fourier transform of the line profile for different wavelength windows. The black line corresponds to the first minimum of the Fourier transform for the wavelength window shown in the right-hand panel. The coloured vertical dashed lines represent the first minimum of the Fourier transform obtained by adjusting the base window by ±0.1 Å the line profile. Right panel: Line profile for the He I+ IIλ4026 line in the stacked spectrum, with the selected wavelength window highlighted in orange.

Determining the precise projected rotational velocity (v sin i) of a rapidly rotating massive star such as ζ Oph can be challenging due to the complexities caused by LPV. Conti & Ebbets (1977) estimated a v sin i of 351 km s⁻¹ using spectra obtained at the Cerro Tololo Observatory and comparing with the non-local thermodynamic equilibrium atmosphere models. Penny (1996) reported a value of 348 km s⁻¹ based on cross-correlation with narrow-lined O-star line-list templates. Howarth & Smith (2001) derived a v sin i of 400 km s⁻¹ using non-local thermodynamic equilibrium atmospheric model fitting, which is consistent with earlier estimates by Herrero et al. (1992) and Balona & Kambe (1999). Jankov et al. (2000) applied the Fourier technique to an averaged spectrum from 242 epochs assembled at the Laboratorio Nacional de Astrofisica, reporting v sin i = 390 km s⁻¹. However, none of these studies accounted for (time-variable) macroturbulence in the broadening estimation of the spectra. The most recent estimate of 319 km s⁻¹ by Simón-Díaz & Herrero (2014) and 385 km s⁻¹ by Holgado et al. (2022) combines the Fourier method with line-profile fitting, explicitly accounting for macroturbulence, but did not consider its time-dependent nature owing to the limited number of epochs.

We determined the v sin i of ζ Oph using the Fourier transform method, following the approach in which the first zero of the Fourier transform of a spectral line profile corresponds to v sin i for stars dominated by rotational broadening (Gray 2005; Simón-Díaz & Herrero 2007, 2014). While generally reliable and a widely used approach, this method can introduce scatter in the measurements, particularly in stars with significant pulsational broadening and/or LPV, which is an important consideration in the spectroscopic analysis of massive stars (see e.g. Aerts et al. 2009, 2014; Simón-Díaz et al. 2017). For stars with significant LPV, this method can lead to uncertainties of up to 5–10% in the inferred v sin i (Piters et al. 1996; Groot et al. 1996). The uncertainty increases if additional broadening, comparable in magnitude to the rotational broadening, is present in the line profile. On the other hand, for fast rotators such as ζ Oph, where rotational broadening dominates and often masks pulsation-induced asymmetries in the line profile, the Fourier method is considered reliable (Simón-Díaz & Herrero 2007; Simón-Díaz et al. 2017).

As an example of the application of the Fourier method to our ESPRESSO data, Fig. 3 illustrates the Fourier transform of the He Iλ4922 line calculated based on the average spectrum of stacking all individual epochs. The line profile is initially selected within a base wavelength window as shown in Fig. 3, then systematically varied by ±0.1 Å at each wing to probe how much our results are sensitive to selection biases and normalisation. The corresponding v sin i values for each iteration are indicated by vertical dashed lines in Fig. 3. To ensure consistency in our subsequent analysis, the line region is manually identified in the stacked spectrum and fixed across all epochs before applying the Fourier transform to derive v sin i. The bootstrapped uncertainty on v sin i measurement from individual epochs is systematically less than 2 km s⁻¹. However, as discussed previously, an uncertainty of ≤5 km s⁻¹ (i.e. ∼3%) is likely underestimated for a star with significant LPV such as ζ Oph (see Piters et al. 1996). The peak-to-peak scatter in the v sin i time series is shown in Fig. 4 for the same three helium lines and is of order 15 km s⁻¹ (i.e. 5% of the mean v sin i value). This scatter can be treated as the lower limit for the uncertainty on v sin i from the Fourier method for a massive star with LPV such as ζ Oph. This is in agreement with the findings of Piters et al. (1996), who suggested that the Fourier method for stars with LPV can yield up to 10% uncertainty of the inferred v sin i. Hence, 15 − 30 km s⁻¹ (i.e. 5–10%) is likely a more realistic uncertainty for our v sin i measurements.

Fig. 4. Left panel: v sin i time series for different spectral lines with peak-to-peak scatter up to 15 km s−1. Right panel: Kernel density estimation of the scatter used to estimate the formal uncertainty in v sin i.

Figure 4 shows the v sin i time series for three different spectral lines, He I+ IIλ4026, He IIλ4200, and He Iλ4922, yielding mean values of 324 ± 13 km s⁻¹, 334 ± 9 km s⁻¹, and 322 ± 13 km s⁻¹, respectively. The indicative formal errors on v sin i from our Fourier analysis for each line are calculated as the 3-σ confidence intervals of the kernel density based on all individual measurements (see Fig. 4), thus accounting for the epoch-to-epoch scatter due to LPV. Again, we stress that these errors might be underestimated, and an uncertainty at a 10% level might be a more reasonable error estimate for stars with LPV following the arguments described above. Interestingly, we find that the v sin i measured from the He IIλ4200 line is systematically larger by about 10 km s⁻¹ compared to the other two lines. The difference in mean v sin i values across the three helium lines could be related to the limitations of the Fourier method as discussed previously. We note that further investigation is necessary to determine the cause of this discrepancy.

We applied our iterative pre-whitening method to the v sin i time series of ζ Oph and detected the 3.33 h, 2.43 h, and 4.63 h periods amongst all three helium lines. As previously mentioned, the 2.43 h and 3.33 h periods are known non-radial pulsation modes, while the 4.63 h period is suggested to be either a radial first overtone or a low-angular degree non-radial mode (Kambe et al. 1990, 1997; Balona & Kambe 1999; Walker et al. 2005). Additionally, all three helium lines exhibited significant periods of 2.91 h and 3.66 h, which were not known in the literature. Furthermore, a period of 7.28 h was found in both the He I + IIλ4026 and He Iλ4922 lines, while the He I + IIλ4026 and He IIλ4200 lines showed periods of 5.37 h and 5.15 h, respectively. These periods are fairly close to the previously detected 5.34 h period, which was suggested to be a low-angular degree non-radial pulsation period by Kambe et al. (1997). Therefore, our analysis demonstrates not only that pulsation periods can be detected in the v sin i time series using the Fourier method, but that the scatter in v sin i for different spectral lines depends on the pulsation properties of the star. Finally, our results highlight the sensitivity of the Fourier method to measuring v sin i in the presence of LPV caused by pulsations in massive stars.

3.3. Macroturbulence

Fig. 5. Macroturbulent velocity (vmacro) for three helium lines in ζ Oph, with the x-axis normalised to the time stamp of the first spectrum (i.e. MJD = 60468.983771). Top panel: time series of all 1706 spectra with 30 s exposure times (in grey). A peak-to-peak variability in vmacro of up to 88 km s−1 is seen. Dashed blue lines are the mean of the vmacro values for each night, which are consistent within their uncertainties as indicated by the shaded region. Bottom panel: vmacro time series for stacked spectra emulating 1 h (black), 30 min (blue), and 15 min (red) integration times. The peak-to-peak variability is diminished for longer integration times, but the mean vmacro is higher for longer integration times, with the nightly mean values denoted as solid lines.

We estimated the macroturbulent velocity broadening (i.e. v_macro) for the three helium lines by assuming a constant rotation rate fixed at 327 km s⁻¹, which is the mean of the v sin i values for all three helium lines, as derived in Sect. 3.2 and shown in Fig. 4. We took the standard approach to constraining macroturbulence in ζ Oph, following the method of Simón-Díaz & Herrero (2007), which assumes that the intrinsic line profile is Gaussian, a full width at half maximum determined by the spectral resolving power, and a line depth set by the observed equivalent width of the line profile. The overall spectral line profile is fitted by including two broadening components – rotation and v_macro – using the isotropic radial-tangential model from Gray (2005). In this isotropic model of macroturbulence, equal fractions of the stellar surface contribute to radial and tangential motion. Whilst this may not be physical, it is a common assumption for massive stars (Ryans et al. 2002; Simón-Díaz et al. 2010, 2017). Following the formalism of Gray (2005), the macroturbulent broadening profile is given by

$$\begin{array}{c}\hfill M(\mathrm{\Delta }\lambda )=\frac{2A_{\mathrm{R}}\mathrm{\Delta }\lambda }{\sqrt{\pi }{\zeta }_{\mathrm{R}}^2}{\displaystyle {\int }_0^{{\zeta }_{\mathrm{R}}/\mathrm{\Delta }\lambda }}{e}^{-1/u^2}\mathrm{d}u+\frac{2A_{\mathrm{T}}\mathrm{\Delta }\lambda }{\sqrt{\pi }{\zeta }_{\mathrm{T}}^2}{\displaystyle {\int }_0^{{\zeta }_{\mathrm{T}}/\mathrm{\Delta }\lambda }}{e}^{-1/u^2}\mathrm{d}u,\end{array}$$(2)

where ζ_R and ζ_T represent the radial and tangential velocity amplitudes, and A_R and A_T are the fractional flux contributions from the radial and tangential components, respectively.

We intentionally adopted an isotropic radial-tangential model for v_macro for our analysis of ζ Oph, which assumes equal flux contributions (i.e. A_R = A_T = 0.5) and equal amplitudes for the radial and tangential velocity components (i.e. ζ_R = ζ_T = v_macro), since this is commonly implemented in spectroscopic studies of massive stars. While these assumptions may not strictly hold, disentangling both the flux ratio and the velocity components from spectroscopy alone is inherently degenerate and necessitates simplifications (e.g. Serebriakova et al. 2023). Moreover, the primary aim here is to demonstrate how common approaches of inferring macroturbulent velocities are significantly affected by LPV caused by pulsations on short timescales. For simplicity, we also fixed the v sin i as mentioned previously, and do not consider the effect of limb darkening or gravity darkening.

We applied this isotropic radial-tangential model to each individual ESPRESSO spectrum of ζ Oph with a 30-s exposure time. Figure 5 shows the v_macro time series for the three helium lines considered in this work. The v_macro measurements from individual 30-s exposures show significant peak-to-peak scatter, up to 88 km s⁻¹, which is much larger than the variability seen in RV or v sin i (cf. Figs. 2 and 4). This clearly indicates that macroturbulence is not a static velocity field, but instead exhibits large-amplitude and short-period variability, which has a specific average value and different temporal behaviour for each spectral line.

We also applied the isotropic radial-tangential model to the time series of stacked spectra that emulate longer exposure times. We did this to investigate potential biases in inferring macroturbulence when integration times are longer and thus more likely to suffer from pulsation-induced smearing of the spectral line profile. Such longer integration times of order tens of minutes are typically needed in high-resolution spectroscopic studies of fainter massive stars. We show the v_macro values derived from stacked spectra emulating 15 min, 30 min, and 1 h integration times in the lower panels of Fig. 5. Importantly, we find that the mean v_macro values for longer integration times are systematically larger, which is due to the pulsational smearing of the line profile for longer exposures, changing the total equivalent width of the line profile.

Moreover, our analysis demonstrates that the peak-to-peak scatter in the v_macro values is line-dependent. For the He IIλ4200 line, the peak-to-peak scatter is systematically larger for shorter integration times, increasing from 40 km s⁻¹ for the 1 h stacked spectra to 88 km s⁻¹ for the individual spectra with 30-s exposure times. For the other two helium lines, the scatter is approximately 20 km s⁻¹ for the 1 h stacked spectra, and exceeds 60 km s⁻¹ for the individual spectra with 30-s exposure times. In Fig. 6 we present the scatter in v_macro as a function of different time-binning, and show that it is systematically larger for the shorter integration times. This implies that high-resolution spectroscopy of fainter massive stars, for which longer integration times are typically required, leads to biased estimates of both the amplitude and variability of macroturbulence. Specifically, the mean macroturbulence may be overestimated due to pulsational smearing, while its intrinsic variability may be underestimated.

Fig. 6. Scatter in the peak-to-peak variability of vmacro for the three helium lines for the different sets of binned spectra.

We applied our iterative pre-whitening method to the v_macro time series for the individual spectra with 30-s exposure times and found that all three helium lines exhibit multi-period variability. The significant periods of 1.39 h, 2.91 h, 3.33 h, 4.63 h, and 5.34 h were observed in all three helium lines. As previously mentioned, the 3.33 h period is a known non-radial pulsation period with spherical harmonic geometry of ℓ = |m| = 4 (Reid et al. 1993; Kambe et al. 1997). The 5.34 h period was also observed by Kambe et al. (1997). Based on the fact that it was only dominant in the line wings and its phase decreased slightly towards longer wavelength, Kambe et al. (1997) suggested that this period is likely due to low-angular degree non-radial pulsation. The 4.63 h period is also likely to be either a radial first overtone or a low-spherical order non-radial pulsation, as reported by Walker et al. (2005).

Additionally, the He I + IIλ4026 and He IIλ4200 lines show periods of 2.43 h and 2.44 h, respectively. These periods are consistent with the previously reported 2.43 h period by Kambe et al. (1993) and Reid et al. (1993). But, whether this period is a true period of a non-radial pulsation mode or an alias of the reported 2.01 h period is debated in the literature kambe1997multiperiodicity; balona1999moving. However, in our analysis, we observed significant periods of 2.04 h and 2.15 h in the He I + IIλ4026 and He IIλ4200 lines, respectively. Within their uncertainty, both these periods are consistent with the 2.01 h period, which is suggested to be a non-radial pulsation with spherical geometry ℓ = 7 − 8 by Kambe et al. (1997) and Balona & Kambe (1999). In summary, our results strongly support that the large-scale and short-period variability of macroturbulence in ζ Oph is caused by multi-periodpulsations.

4. Impact of variable macroturbulence on inferred atmospheric parameters

Having established the significant large-scale and short-period multi-period variability in macroturbulence for ζ Oph, which is much larger than that of the RV and v sin i time series, we also investigated the impact this has on inference of stellar parameters. We used the grid-based quantitative spectroscopy tool IACOB-GBAT (Simón-Díaz et al. 2011), which employs FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Sundqvist & Puls 2018; Puls et al. 2020) atmospheric models to derive T_eff and log g for subsets of spectra that sample the full range in v_macro as shown in Fig. 5. In this way, the full range of v_macro values are propagated to create distributions for T_eff and log g and test by how much these parameters vary because of variable macroturbulence.

Figure 7 shows histograms of the derived T_eff and log g estimates for a subset of 100 spectra with 30-s exposure times. We find a non-negligible scatter in the atmospheric parameters caused by variable macroturbulence, with a range almost as large as the formal uncertainties for each parameter when considering a single spectrum. This demonstrates that while time-dependent macroturbulence contributes to systematic and random scatter in the inferred stellar parameters, such effects are only detectable with the ultra-high precision offered by high-resolution spectrographs such as ESPRESSO for the brightest massive stars.

Fig. 7. Histograms of Teff and log g for the 30 s spectra. The vertical dashed lines in black indicate the best-fit values obtained from the full-night stacked spectra. The red and blue arrows indicate the average confidence interval on the Teff and log g estimates from individual spectra.

5. Comparison to archival and TESS photometry

Fig. 9. Top panel: extracted TESS light curve from sector 91 of ζ Oph extracted using simple aperture photometry. Bottom panel: Lomb-Scargle periodogram of the TESS light curve, with the significant periods identified in this work shown as dashed grey lines.

Fig. 10. From top to bottom: Lomb-Scargle periodograms for the RV, v sin i, and vmacro time series for the same three helium lines. In each panel, the dashed blue lines indicate the significant periods extracted for each spectroscopic time series using the pre-whitening method described in Sect. 3, the black circles indicate the pulsation periods extracted from the TESS light curve, and the red triangles correspond to the pulsation periods from ground-based and MOST photometry reported by Walker et al. (2005).

ζ Oph has been studied extensively using ground-based photometry as well as space-based photometry from the MOST, SMEI, and WIRE missions. Several independent analyses have revealed multi-period pulsations with quasi-periodic amplitudes and frequencies (e.g. Balona & Kambe 1999; Walker et al. 2005; Howarth et al. 2014). All photometric studies consistently report a dominant pulsation period of 4.63 h. Using MOST photometry, Walker et al. (2005) also detected the spectroscopically identified periods of 2.43 h, 3.33 h, and 5.34 h, but did not recover the 2.01 h period reported previously by Kambe et al. (1997). Following the approach of Telting & Schrijvers (1997) and using simultaneous photometry and spectroscopy, Walker et al. (2005) suggested that the 4.63 h pulsation period is a non-radial mode with ℓ ∈ {3, 4}, while the 2.43 h and 3.33 h periods correspond to ℓ = 7 and ℓ ∈ {3, 4}, respectively. However, using non-adiabatic pulsation analysis for evolutionary models of a 23-M_⊙ star, Walker et al. (2005) found that the 4.63 h period may also be consistent with a radial first overtone or a non-radial mode with (ℓ,m) = (1, −1) or (3, 1).

The disagreement in the inferred mode geometry between different studies, or between methods within the same study, may arise because of coupling between rotation and pulsation modes in ζ Oph, which means that the geometry of a pulsation mode cannot be described by a single spherical harmonic. For example, coupling between different spherical harmonics, such as between spheroidal and toroidal components, can occur (see Telting & Schrijvers 1997). This can lead to different inferred ℓ values depending on whether such coupling is considered or not when modelling spectroscopic or photometric variability (see e.g. Daszyńska-Daszkiewicz et al. 2002). Therefore, whilst it is clear that ζ Oph is multi-period pulsator, there is still ambiguity in the spherical harmonic geometries of its pulsation periods.

To complement our spectroscopic analysis of ζ Oph, we also analysed recently assembled TESS space photometry. The NASA TESS mission is an all-sky survey providing optical light curves at high cadence and had a nominal mission to identify Earth-like planets transiting bright stars (Ricker et al. 2015). TESS observed ζ Oph for the first time in its sector 91 in April 2025. Despite a minor safe-mode event, TESS assembled 28-d time-series photometry with a cadence of 200 s in its full-frame images (FFIs). We used the LIGHTKURVE software package (Lightkurve Collaboration 2018) to take a cutout of the TESS FFIs for sector 91 and followed the methodology of Bowman et al. (2022) to select an optimum aperture mask to extract a light curve. We then used principal component analysis to detrend the light curve of remaining instrumental trends. Since ζ Oph is a bright star, it is saturated in the TESS FFIs, but we were able to capture sufficient flux by including the bleed columns within our aperture mask, which is shown in Fig. A.1. The spectral window of the extracted TESS light curve is shown alongside that of the ESPRESSO time series in Fig. 8, which demonstrates the improved nature of aliasing for near-continuous space photometry compared to ground-based spectroscopy.

The detrended sector 91 TESS light curve and its corresponding Lomb-Scargle periodogram are shown in Fig. 9. We applied an iterative pre-whitening method as before to detect the significant periods from the TESS dataset, which are indicated in Fig. 9 and included in Table B.1 along with their least-squares fitting uncertainties. These new TESS data allow us to recover several significant pulsation periods that have been previously reported in ground-based photometric studies by Balona & Kambe (1999) and Walker et al. (2005), but also detect several new significant pulsation periods.

To compare the significant periods found in the TESS light curve and the spectroscopic RV, v sin i, v_macro time series, we show the Lomb-Scargle periodograms of each spectroscopic dataset in Fig. 10. For ease of comparison between the datasets, significant periods from the TESS light curve are marked by black circles, while periods from archival photometry (Walker et al. 2005) are marked with red triangles, and significant periods from each spectroscopic time series are shown as vertical dashed blue lines. Among these, the previously known pulsation periods of 3.33 h and 2.43 h are consistently detected across all spectroscopic datasets and are also clearly present in the TESS light curve (also see Kambe et al. 1993, 1997; Balona & Kambe 1999). A period of 2.91 h is similarly detected in both photometric and spectroscopic data, suggesting that it represents a genuine pulsation frequency, although this period has never been reported before in the literature.

Additionally, the 4.63 h period, which is significant in both the v sin i and v_macro time series, is also identified in the TESS photometry with the highest S/N of all pulsation periods. This periodicity was initially reported by Balona & Kambe (1999) based on ground-based observations from 1985 but was notably absent in subsequent photometric campaigns. However, Walker et al. (2005) and Howarth et al. (2014) again reported the detection of the 4.63 h period based on the MOST, SMEI and WIRE photometry from 2003 to 2008. Its detection in the TESS light curve with high significance strongly supports the interpretation that the 4.63 h period corresponds to a pulsation period rather than an alias.

The 5.34 h pulsation period reported by Kambe et al. (1997) and Walker et al. (2005) is also recovered in the TESS data. We also detect significant periods of 5.42 h, 5.43 h, and 5.36 h in the v_macro time series of the three diagnostic helium lines from our ESPRESSO spectra, all of which are consistent with the 5.34 h period within their uncertainties. The 2.01 h pulsation period is not detected in the TESS photometry, despite being consistently observed in our RV time series for all three helium lines, as well as in the v sin i and v_macro time series for the He I +IIλ4026 line and the v sin i time series of the He Iλ4922 line. Notably, the sum of the frequencies corresponding to the periods of 4.63 h (f₁) and 3.56 h (f₂) (i.e. a combination frequency of f₁ + f₂) yields a frequency at 11.93 d⁻¹ in the TESS data, which is very close to the 11.98 d⁻¹ frequency associated with the reported alias period of 2.01 h. This supports the interpretation that the 2.01 h period is a combination frequency, while the 2.43 h period is the real pulsation period rather than an alias of the 2.01 h period.

While several other longer periods known from archival photometry were also observed in the TESS data, these were not observed in the ESPRESSO data. However, this is because the ESPRESSO spectroscopic time series has a significantly different temporal sampling and spectral window, resulting in reduced sensitivity at low frequencies. For example, the time span of our ESPRESSO data is four consecutive nights, meaning that periods longer than several days are subject to aliasing. Consequently, these low-frequency signals cannot be identified with confidence in the spectroscopic diagnostics.

6. Discussion and conclusions

Using high-cadence, high-resolution, and high-S/N ESPRESSO spectroscopy, we provide compelling evidence that macroturbulence in the rapidly rotating and pulsating massive star ζ Oph is not a static velocity field but instead exhibits large-amplitude, short-period, and line-specific variability caused by pulsations. The v_macro time series derived from individual 30 s spectra in particular shows a peak-to-peak scatter of up to 88 km s⁻¹ and periods of order hours. Furthermore, using different stacked time series to emulate longer integration times, we demonstrate that the inferred amplitude of the variability in macroturbulence depends on the integration time of the spectra because of a pulsation-induced smearing of the line profile. For spectra with integration times comparable to the shortest pulsation period (i.e. ∼1 h), the observed scatter in v_macro is significantly reduced down to ∼20 km s⁻¹. Moreover, the mean v_macro values inferred from longer integration times are systematically higher than those derived from shorter integration time spectra. This is due to the pulsational smearing of the line profile in long-exposure observations. This result is similar to previous works by Kambe et al. (1997) and Aerts et al. (2009), who looked at the artificial broadening of spectral lines caused by pulsations in massive stars. This suggests that traditional approaches that rely on a single stacked spectrum systematically overestimate macroturbulent broadening.

However, there are limitations to our analysis. The helium lines we used do have small amounts of Stark broadening, and hence our estimations of macroturbulent velocities are slightly overestimated. However, the Stark broadening is much smaller than the total non-rotational additional broadening we refer to as macroturbulence. Additionally, the assumption of isotropic macroturbulence does not fully hold when the observational cadence is much shorter than the dominant pulsation period. In such cases, pulsations can produce asymmetric line profiles, and the symmetric isotropic macroturbulence may not be an accurate representation of the line profile. This effect is relatively small compared to the dominant rotational broadening for rapidly rotating stars such as ζ Oph, but remains non-negligible. Regardless, our study has concretely demonstrated how conventional methods of inferring macroturbulence can be biased by pulsations and instrumental setup. As future work, we will model the pulsational effects in more detail by incorporating line asymmetry and accounting for the differing contributions of radial and tangential velocity components from pulsations. Disentangling the ratio of radial to tangential velocities is particularly important, as gravity modes typically have ζ_t ≫ ζ_r, which may differ from other instabilities such as convection.

We assessed the impact of time-dependent macroturbulence on stellar parameter inference by applying the grid-based fitting tool IACOB-GBAT to spectra, sampling the full range of observed v_macro values. This approach enabled us to construct distributions for T_eff and log g. The scatter in the derived parameters is comparable to the formal uncertainties on an individual spectrum and indicates that macroturbulence variability can introduce both systematic and random errors in parameter estimation. However, these effects are only apparent with the ultra-high S/N and time resolution of ESPRESSO and would likely remain undetected in conventional, time-averaged datasets.

Iterative pre-whitening applied to the RV, v sin i, and v_macro time series revealed a large number of significant periods, further evidence that ζ Oph is a multi-period pulsator with significant LPV. Using the first TESS light curve ever assembled for ζ Oph, we also identified the significant pulsation periods, many of which are also found in our ESPRESSO spectroscopic time series. The 2.43 h, 2.91 h, 3.33 h, and 4.63 periods were observed consistently in both the v sin i and v_macro time series across all three diagnostic helium lines, and in our TESS light curve. Previous spectroscopic investigations by Kambe et al. (1990, 1993, 1997) and Balona & Kambe (1999), along with photometric follow-up studies by Walker et al. (2005) and Howarth et al. (2014), provided strong evidence that the 3.33 h and 2.43 h periods are caused by non-radial pulsations. A period of 5.34 h was found in both the TESS light curve and the v_macro time series; this was also found in previous studies and suggested to be a low-angular-degree non-radial pulsation (Kambe et al. 1997). However, there remains some debate in the literature as to which of the 2.01 h or 2.43 h periods is an alias. In our analysis, the 2.43 h period is clearly detected in the TESS photometry, while the 2.01 h signal is not; thus, we conclude that, given TESS’s excellent duty cycle and spectral window, the 2.01 h period is the alias.

The similarity of the significant periods in our spectroscopic data and those identified as pulsation modes in TESS and archival photometry provides compelling evidence that the observed variability in macroturbulence in ζ Oph is driven by non-radial pulsations. Thus, our proof of concept using ESPRESSO data of ζ Oph demonstrates that pulsations play a critical role in the physics of (time-dependent) macroturbulence in massive stars. The variability observed in v_macro is not only large in amplitude but also quasi-periodic and spectral-line-specific, and we demonstrate how it can be overestimated for fainter massive stars that require long exposure times. However, there remains ambiguity in the identification of spherical harmonic geometry for ζ Oph’s pulsation modes, which will be the subject of our future work. Incorporating significant frequencies identified from TESS and archival photometry into physically motivated, time-dependent, and anisotropic prescriptions for macroturbulence will be instrumental in further improving our understanding of the structures of massive stars.

Acknowledgments

The authors gratefully acknowledge UK Research and Innovation (UKRI) in the form of a Frontier Research grant under the UK government’s ERC Horizon Europe funding guarantee (SYMPHONY; PI Bowman; grant number: EP/Y031059/1), and a Royal Society University Research Fellowship (PI Bowman; grant number: URF/R1/231631). This project received the support from the “La Caixa” Foundation (ID 100010434) under the fellowship code LCF/BQ/PI23/11970035. SS-D acknowledges support from the State Research Agency (AEI) of the Spanish Ministry of Science and Innovation (MICIN) and the European Regional Development Fund, FEDER under grants LOS MÚLTIPLES CANALES DE EVOLUCIÓN TEMPRANA DE LAS ESTRELLAS MASIVAS PID2021-122397NB-C2. Based on observations collected with the ESPRESSO spectrograph at the European Southern Observatory in Paranal, Chile, under ESO programme 113.26B9 (PI Bowman), which are accessible through the ESO archive (https://www.eso.org/sso/login). Data products that support the results in this paper are publicly available via the Zenodo repository: https://doi.org/10.5281/zenodo.17233369.

References

Appendix A: TESS photometry

In Fig. A.1 we provide the TESS aperture and background masks used to extract the light curve in Fig. 9.

Fig. A.1. Our custom aperture and background masks based on the TESS FFI image of ζ Oph in sector 91.

Appendix B: Extracted periodicity

Significant periods identified across different time series, including from TESS and archival photometry, and ESPRESSO spectroscopy, and their respective uncertainties are provided in Table B.1.

All Tables

All Figures

	Fig. 1. Left panels: three spectra from 9 June 2024 showing LPV for the He I+IIλ4026, He IIλ4200, and He Iλ4922 lines (from top to bottom), with the stacked mean spectrum overplotted in dark blue. Right panels: residual dynamic spectra for the corresponding left panels.
In the text

	Fig. 2. RV time series for three helium lines in ζ Oph extracted using cross-correlation. The peak-to-peak RV variations exceed 20 km s−1 in a multi-period manner because of pulsations rather than binarity.
In the text

	Fig. 8. Spectral windows for the ESPRESSO spectroscopic time series (top panel) and the TESS photometric data (bottom panel) analysed in this work.
In the text

Fig. 3. Left panel: Fourier transform of the line profile for different wavelength windows. The black line corresponds to the first minimum of the Fourier transform for the wavelength window shown in the right-hand panel. The coloured vertical dashed lines represent the first minimum of the Fourier transform obtained by adjusting the base window by ±0.1 Å the line profile. Right panel: Line profile for the He I+ IIλ4026 line in the stacked spectrum, with the selected wavelength window highlighted in orange.

In the text

	Fig. 4. Left panel: v sin i time series for different spectral lines with peak-to-peak scatter up to 15 km s−1. Right panel: Kernel density estimation of the scatter used to estimate the formal uncertainty in v sin i.
In the text

Fig. 5. Macroturbulent velocity (vmacro) for three helium lines in ζ Oph, with the x-axis normalised to the time stamp of the first spectrum (i.e. MJD = 60468.983771). Top panel: time series of all 1706 spectra with 30 s exposure times (in grey). A peak-to-peak variability in vmacro of up to 88 km s−1 is seen. Dashed blue lines are the mean of the vmacro values for each night, which are consistent within their uncertainties as indicated by the shaded region. Bottom panel: vmacro time series for stacked spectra emulating 1 h (black), 30 min (blue), and 15 min (red) integration times. The peak-to-peak variability is diminished for longer integration times, but the mean vmacro is higher for longer integration times, with the nightly mean values denoted as solid lines.

In the text

	Fig. 6. Scatter in the peak-to-peak variability of vmacro for the three helium lines for the different sets of binned spectra.
In the text

	Fig. 7. Histograms of Teff and log g for the 30 s spectra. The vertical dashed lines in black indicate the best-fit values obtained from the full-night stacked spectra. The red and blue arrows indicate the average confidence interval on the Teff and log g estimates from individual spectra.
In the text

	Fig. 9. Top panel: extracted TESS light curve from sector 91 of ζ Oph extracted using simple aperture photometry. Bottom panel: Lomb-Scargle periodogram of the TESS light curve, with the significant periods identified in this work shown as dashed grey lines.
In the text

Fig. 10. From top to bottom: Lomb-Scargle periodograms for the RV, v sin i, and vmacro time series for the same three helium lines. In each panel, the dashed blue lines indicate the significant periods extracted for each spectroscopic time series using the pre-whitening method described in Sect. 3, the black circles indicate the pulsation periods extracted from the TESS light curve, and the red triangles correspond to the pulsation periods from ground-based and MOST photometry reported by Walker et al. (2005).

In the text

	Fig. A.1. Our custom aperture and background masks based on the TESS FFI image of ζ Oph in sector 91.
In the text