Possible anti-correlations between pulsation amplitudes
and the disk growth of Be stars in giant-outbursting Be X-ray binaries

Neil Gehrels Swift Observatory (Swift; Gehrels et al., 2004) is a gamma-ray burst survey satellite operated by NASA. We downloaded the Burst Alert Telescope (BAT) 15–50 keV daily light curves of target BeXBs, from Swift/BAT Hard X-ray Transient Monitor web-page (Krimm et al., 2013). Seven days binning was performed to improve signal-to-noise ratio for each light curves.

Monitor of All-sky X-ray Image (MAXI; Matsuoka et al., 2009) is an X-ray camera installed onboard the International Space Station (ISS). We obtained weekly-averaged Gas Slit Camera (GSC) 2–20 keV light curves covering 58000–60000 MJD for target BeXBs, via MAXI on-demand web interface (Nakahira et al., 2012).

Zwicky Transient Facility (ZTF; Bellm et al., 2019) is an optical wide-field survey project using a CCD camera array attached to Samuel Oschin telescope at Palomar Observatory in California, operated by California Institute of Technology (Caltech). We utilised $g^{\prime}$ and $r^{\prime}$-band light curves of the targets acquired from ZTF Public Data Release 15 (DR15; Masci et al., 2019), which covers the epoch of $\sim$ 58250–59900 MJD. The light curve is not available only for RX J0440.9+4431 due to its brighter optical magnitude than the ZTF saturation magnitude ($\sim$ 12.5 AB).

Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al., 2018) is a high-cadence optical all-sky survey project using four 50 cm telescope systems located in Haleakala, Mauna Loa (Hawaii), El Sauce (Chile), and Sutherland (South Africa). We obtained $c$ and $o$-band light curves covering 58000–60000 MJD via ATLAS Forced Photometry web interface (Smith et al., 2020). The saturation magnitude of ALTAS is $\sim$ 12.5 AB, about same as ZTF. Thus, the light curve for RX J0440.9+4431, although it can be created, is not reliable.

Gattini-IR (De et al., 2020) is a NIR wide-field survey project using 30 cm robotic $J$-band telescope system at Palomar Observatory, operated by Caltech. Gattini-IR sweeps $\sim$ 7500 $\mathrm{deg}^{2}$ of the sky at one night with an upper-limit magnitude of $\sim$ 15.7 (AB), and scans the northern-sky observable from Palomar (Dec. $\gtrsim-30\degr$) with a cadence of $\sim$ 2 days. The system is designed for the NIR-transient detection, and a real-time data processing pipeline (Gattini Data Processing System; GDPS) analyses observed images and detects the transients. We obtained the $J$-band light curves of the targets generated by GDPS. The data-set covers the epoch of $\sim$ 58400–59900 MJD.

The data reduction process consists of three stages: detrending, 5-sigma clipping, and 30-min binning. First, we subtracted the moving averaged flux from the light curve to remove the long-term flux variation in time scale of $\gtrsim$ several days. For convenience, we call this process as ‘detrending’. The width of the detrending window was set to 3 days except for SAX J2103.5+4545. For SAX J2103.5+4545, the width was set to 1 day because the 3-day detrending was insufficient to remove long-term variability. Next, 5-sigma clipping was performed to remove outliers. Finally, we binned light curves in 1/48 days ($\sim$ 30 min) width to minimise the effect of differences in the sampling rate of FFIs between the prime and extended missions (Ricker, 2021). Figure 3 shows the data reduction process to the TESS light curve of V0332+53, sector 59.

We applied the Lomb-Scargle technique to the reduced light curves. Note that we did not use a power spectrum (periodogram), but an amplitude spectrum calculated by the following transformation (cf. Aerts, 2021):

where $\nu$ is a frequency, $A(\nu)$ is the amplitude spectrum, $N$ is a number of points in the light curve, and $P(\nu)$ is the power spectrum. We first made dynamic amplitude spectra. A dynamic amplitude spectrum is obtained by taking a portion of the light curve with an appropriate window and obtaining its spectrum by applying the Lomb-Scargle while moving the window. That shows the variation of the spectrum. The window width was set to 5 days and moved by 1/48 days. Note that oscillations with a period equal to or greater than the width of the detrending window had removed in the data reduction process. Figure 4 is a dynamic amplitude spectrum of V0332+53, sector 59.

If the number of points in the extracted light curve was less than 160, the obtained spectrum was considered to be unreliable, where the number 160 is 2/3 of the number of points expected from the sampling rate (48 $\mathrm{d^{-1}}$) and the window width (5 days). At the same time, we evaluated the contribution of photometric errors to the spectrum by creating mock light curves which fluctuate randomly with a Gaussian-distribution over the width of the photometric error and calculating their spectrum. Note that the photometry pipeline we used has a tendency to underestimate errors by about 10 per cent and up to 50 per cent, compared with estimated from the light curve fluctuations. Finally, we averaged each dynamic spectrum along the time-axis direcion (Figure 5).

In this calculation, spectrum with points fewer than 160 on the light curve were excluded. We defined this as "averaged amplitude spectrum", and considered it to be the representative amplitude spectrum of the object in the sector. The reason why Lomb-Scargle is not simply applied to the entire light curve for each sector is to reduce the effect of several days of missing points on some light curves.

To extract the peaks from the spectrum, we estimated a red-noise component of the spectrum. First, we made a local 25th percentile (first quartile) curve by extracting a portion of the spectrum with a window moving across all frequencies, and calculating the 25th percentile for each window. Then the 25th percentile curve is smoothed with a Gaussian filter. The obtained curve can be considered to be a rough estimate of the continuum. The reason for using the 25th percentile rather than the mean or median is the nature of the amplitude spectrum. Unlike random errors, which can be both positive and negative, the peaks in the amplitude spectrum are almost always positive. Therefore, the mean or median are expected to be systematically larger than the continuum. Finally, the red noise spectrum was obtained by model fitting to the roughly estimated continuum. We utilised the model used in Bowman et al. (2019) and Nazé et al. (2020) with some modifications. The model consists of a Lorentzian-like term and a constant term:

where $\nu$ is a frequency, $A_{\mathrm{RN}}$ is a red noise amplitude spectrum, $a_{0}$ is a scaling factor, $\nu_{0}$ is a characteristic frequency, $\gamma$ is an exponent, and $C$ is a white noise level. When $\gamma=2$, the first term is the Lorentzian. If both sides of Equation 8 are squared, the equation agrees with Bowman et al. (2019) and Nazé et al. (2020). The square root is applied to make the dimension of both sides amplitude rather than power. Figure 6 shows the red noise estimation of V0332+53, sector 59.

Incidentally, the error-derived amplitude (cf. Section 4.3.2) was smaller than the white noise level for all light curves as shown in Figure 6, even if putting into account the underestimation of the error.

The peaks were extracted by subtracting the red noise from the amplitude spectrum. We determined that the peaks whose amplitude were larger than twice the red noise were significant. For each source, we evaluated the amplitudes for all sectors at frequencies for which a significant peak was detected in at least one of the sectors. The amplitude of the detected peak was normalised by the flux calculated from the Tmag of the source. This normalisation is based on the assumption that Tmag is the typical magnitude of the source in the TESS pass-band, and ignores long-term optical variations.

In the first scenario, we consider that changes in the density distribution of co-rotating gas clouds are the cause of amplitude variations. In other words, the origin of flux oscillations are not pulsations but the rotation. Since rotation-derived flux oscillations are caused by co-rotating gas clouds, the amplitude of flux oscillation should be determined by the dynamic range of gas cloud density around the rotation axis. Therefore, the peak amplitude is expected to decrease in conjunction with the mass ejection of the Be star and the disk growth due to that the amount of co-rotating gas increased and homogenised around the rotation axis (Figure 23).

However, unlike pulsation, which can have multiple modes, rotation is a single mode oscillation, so it is unlikely that rotation-derived anti-correlated and non-correlated peaks coexist. Furthermore, this scenario provides no explanation for amplitude variations of pulsation-derived flux oscillations. Therefore, in this scenario, all anti-correlated and non-correlated peaks must be rotation-derived and pulsation-derived, respectively. Nevertheless, only 4U0115+634 and RX J0440.9+4431 have peaks with features suggesting a rotation origin, and to apply this scenario to other sources, it would have to be that the harmonics were detected but the fundamental mode was not.

The second scenario is based on the idea that the amplitude of the pulsation itself are constant, but the flux oscillation amplitude varies linked to the change in the amount of circumstellar materials (e.g. disk, stellar wind) and their reprocessing. We considers a situation in which most or all of the photosphere is covered by the circumstellar materials, as seen from our perspective. In this scenario, the photons coming from the Be star are scattered by the circumstellar materials, with at least 20 per cent of the photons passing through without being scattered even once (Figure 24).

If the main scattering process is Thomson scattering, the reprocessing itself does not affect on the SED. However, because photons emitted from various regions on the photosphere are mixed, the scattered flux is considered to have lost periodicity. Since the mass ejected from the Be star should be supplied not only to the disk but also to the circumstellar materials, we expect to observe the peak amplitude decrease in conjunction with the disk growth. This logic can also be applied to explain rotation-derived oscillations. However, all peak amplitudes should decrease uniformly because almost all photons are reprocessed at the same rate in this mechanism. Thus, this scenario cannot explain non-correlated peaks.

Because the situation in this scenario is simple, it is easy to estimate the density of the circumstellar materials. The following relationship is expected to be satisfied when photons pass through a material:

where $\tau$ is a transmittance, $N$ is a electron column density, and $\sigma$ is a scattering cross section. Substituting $\tau=0.2$ and the Thomson scattering cross section $\sigma_{T}=6.7\times 10^{-25}~{}\mathrm{cm^{2}}$ for $\sigma$, the column density is calculated as $N=2.4\times 10^{24}~{}\mathrm{cm^{-2}}$. Assuming a spherical symmetric stellar wind as scattering materials, an electron number density $n$ is expected to follow an inverse square law with respect to a distance from the star $r$ due to the conservation of particle number, when ignoring particle acceleration:

where $n_{0}$ is a number density at the surface and ${R_{*}}$ is a radius of the star. $N$ is then equal to integral of $n$ from the stellar surface to infinity:

If ${R_{*}}=10{R_{\sun}}$, then $n_{0}=3.4\times 10^{12}~{}\mathrm{cm^{-3}}$. Assuming that the material is hydrogen, the numbers of electrons and nucleons are equal, thus a density at the surface is $\rho=6\times 10^{-12}~{}\mathrm{g\,cm^{-3}}$. This is comparable to the observed density of the inner disk ($=10^{-13}\text{--}10^{-10}~{}\mathrm{g\,cm^{-3}}$, cf. Gies et al., 2007; Pott et al., 2010; Schaefer et al., 2010; Štefl et al., 2012; Meilland et al., 2012). Figure 25 shows the estimated density profile.

The estimated density is $\sim$ 100 times larger than the numerical calculated wind density profile of the B1Ve star by Curé (2004). If the density distribution is not spherically symmetric, i.e., if the circumstellar material in our line of sight is locally dense, the less average density can be accounted for. In this case, there should be a clump at least 100 times denser than normal spherical symmetric winds. For example, the winds of Wolf-Rayet (WR) stars are clumpy as indicated by the stochastic variation in the line profiles (Moffat et al., 1988). The flux-limiting radiation hydrodynamic simulation of Moens et al. (2022) predicted a dynamic range of $\sim 100$ in density of clumpy winds of WR stars. In addition, analysis of the UV and optical spectra of eight galactic O-type supergiants by Hawcroft et al. (2021) showed that half of the wind velocity field was covered by dense clumps and that their density was 25 times the average density.

Incidentally, it is not obvious whether the effect of the temperature increase of the photosphere due to the return of some of the scattered photons to the Be star can be ignored in the situation where 80 per cent of the photons are scattered. So, we calculated the fraction of returning photons and obtained the value of about 10–20 per cent in the case considered here (Appendix B). Thus, at least the effect on the estimation of Be star properties in Section 4.1 can be considered negligible.

The third scenario shares the basic idea with the second one (Section 6.2.2) that the circumstellar materials obscure the Be star, but differs in that they obscure only a portion of the star, not the entire star (Figure 26).

In this case, since some areas on the photosphere are not covered, the presence of non-correlated peaks can be tolerated. The question here is whether it is possible to reduce the peak amplitude by 80 per cent just by hiding part of the photosphere. For example, assuming the obscuration by the disk as shown in Figure 26, the obscured area is at most 50 per cent of the entire photosphere, and it seems difficult to achieve the amplitude reduction of 80 per cent. It can be satisfied in the case where the flux oscillations are generated in a limited region around the Be star. For the rotation-derived flux oscillations, this logic is applicable because they originate from co-rotating gas clouds. However, the probability that the co-rotating gas cloud is hidden by the disk is at most 50 per cent, so it is slightly difficult to explain why significant anti-correlations were found in four of the five sources. On the other hand, it is non-trivial whether the requirement is satisfied for the pulsation-derived flux oscillations. Non-radial pulsations cause the local expansion and contraction, and the geometry of the oscillation is characterised by the spherical harmonics. The local brightness at the expanding region is expected to be smaller due to gravitational darkening. In low-$l$ modes, the number of nodal surfaces appearing on the photosphere is less, so each area of expansion or contraction is relatively large, and it is thought that the expanding and contracting regions do not appear equal when viewed from one direction. In this case, periodic flux oscillations can be generated by the brightness variations mentioned above, and the area contributing to the flux oscillations is relatively large. On the other hand, the same logic cannot simply be applied to cases where the expanding and contracting regions appear equal, such as in high-$l$ modes. In this case, flux oscillations can be generated by a local brightness difference (e.g. star spots) between the expanding and contracting regions. Therefore, some of the pulsation-derived flux oscillations may be caused by a limited region on the photosphere. However, the amplitude of the flux oscillations produced by this mechanism is considered to be smaller than that produced by the expansion and contraction of a larger region in low-$l$ modes. That makes difficult to interpret the anti-correlations of the maximum peaks in this scenario. In addition, local brightness differences on the photosphere are expected to be visible and hidden as the Be star rotates, which would cause rotation-linked features in light curves. Therefore, it is also difficult to apply this scenario to SAX J2103.5+4545 and V0332+53, which have no significant rotation-derived features in their amplitude spectra.

In the fourth scenario, we considered that the pulsation amplitude is indeed inversely correlated with the disk growth. The pulsations in early-type massive stars are driven by the $\kappa$-mechanism (Dziembowski & Pamiatnykh, 1993; Dziembowski et al., 1993; Gautschy & Saio, 1993; Pamyatnykh, 1999; Miglio et al., 2007). In the temperature range where the Rosseland opacity $\kappa$ is positively correlated with the temperature $T$, $\kappa$ variations provide positive feedback to $T$ variations, causing growth of thermodynamic fluctuations and thus pulsations. Since the amount of this feedback should depend on the slope of the $T\text{-}\kappa$ curve (i.e. $\mathrm{d}\kappa/\mathrm{d}T$), the intensity of $\kappa$-mechanism, and thus the pulsation amplitude, is expected to have dependence on the inner temperature distribution. Therefore, we first assume that some variation in the temperature inside the Be star has caused the pulsation to become inactive. The $\kappa$-mechanism is expected to be most active at the steepest point on the low-temperature side slope of the iron-group element derived bump (so-called Z-bump) in the $T\text{-}\kappa$ curve. And therefore, $\mathrm{d}\kappa/\mathrm{d}T$ becomes smaller in the higher temperature as it is closer to the Z-bump peak. Conversely, moving to the low temperature side also reduces $\mathrm{d}\kappa/\mathrm{d}T$ as it approaches the bottom of valley between Z-bump and the hydrogen-derived bump. In other words, a decrease in pulsation amplitude is qualitatively expected for both temperature increases and decreases (Figure 27).

If the star is hotter in the pulsation-inactive state than the pulsation-active state, relatively strong radiation pressure may promote mass ejection. On the other hand, considering the low temperature pulsation-inactive state, the lower temperature may be caused by the star expansion and the outer layer is more loosely bound than the active state, which is also expected to be relatively favourable for mass ejection. In any case, the pulsation-inactive state is qualitatively expected to be more active in mass ejection than the pulsation-active state. Note that these are only qualitative discussions, and a more quantitative discussion requires a detailed understanding of the structure and temperature distribution inside Be stars.

Since this mechanism works for all pulsations, it cannot account for pulsation-derived non-correlated peaks. In addition, this scenario provides no explanation for the rotation-derived flux oscillations. Therefore, the origin of the anti-correlated and non-correlated peaks should be pulsation and the rotation, respectively. This is not so unreasonable since most of the anti-correlated peaks did not have the rotation-derived features, as discussed in Section 6.1.