The mechanism of Xray outbursts in Be Xray binaries remains a mystery and understanding their circumstellar discs is crucial for a solution of the masstransfer problem. In particular it is important to identify the Be star activities e.g. pulsations that cause mass ejection and hence disc formation. Therefore we investigated the relationship between optical flux oscillations and the infrared IR excess in a sample of five Be Xray binaries. Applying the LombScargle technique to highcadence optical light curves from the Transiting Exoplanet Survey Satellite inlineformulatexmath idTM0001 notationLaTeXit TESStexmathinlineformula we detected several significant oscillation modes in the 324 h period range for each source. We also measured the IR excess a proxy for disc growth of those five sources using Jband light curves from Palomar GattiniIR. In four of the five sources we found anticorrelations between the IR excess and the amplitude of the main flux oscillation modes. This result is inconsistent with the conventional idea that nonradial pulsations drive mass ejections. We propose an alternative scenario where internal temperature variations in the Be star cause transitions between pulsationactive and massejectionactive states.
['Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Possible anticorrelations between pulsation amplitudes and the disc growth of Be stars in giantoutbursting Be Xray binaries
{'Possible anti-correlations between pulsation amplitudes and the disk growth of Be stars in giant-outbursting Be X-ray binaries': "Masafumi Niwano, 1 ★ Michael M. Fausnaugh, 2 Ryan M. Lau, 3 Kishalay De, 4 Roberto Soria, 5 , 6 , 7 George R. Ricker, 4 Roland Vanderspek, 4 Michael C. B. Ashley, 8 Nicholas Earley, 9 Matthew J. Hankins, 10 Mansi M. Kasliwal, 9 Anna M. Moore, 11 Jamie Soon, 11 Tony Travouillon, 11 Mahito Sasada, 12 Ichiro Takahashi, 1 Yoichi Yatsu, 1 and Nobuyuki Kawai 1  \n1 Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan  \n2 Department of Physics & Astronomy, Texas Tech University, Lubbock, TX 79410-1051, USA  \n3 NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA  \n- 4 MIT-Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  \n5  \n6  \nCollege of Astronomy and Space Sciences, University of the Chinese Academy of Sciences, Beĳing 100049, China  \nINAF-Osservatorio Astrofisico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy  \n- 7 Sydney Institute for Astronomy, School of Physics A28, The University of Sydney, Sydney, NSW 2006, Australia  \n8  \nSchool of Physics, University of New South Wales, Sydney NSW 2052, Australia  \n- 9 Cahill Center for Astrophysics, California Institute of Technology, 1200 E. California Blvd. Pasadena, CA 91125, USA 10\n- Arkansas Tech University, Russellville, AR 72801, USA\n- 11 Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia\n- 12 Institute of Innovative Research, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'The mechanism of X-ray outbursts in Be X-ray binaries remains a mystery, and understanding their circumstellar disks is crucial for a solution of the mass-transfer problem. In particular, it is important to identify the Be star activities (e.g., pulsations) that cause mass ejection and, hence, disk formation. Therefore, we investigated the relationship between optical flux oscillations and the infrared (IR) excess in a sample of five Be X-ray binaries. Applying the Lomb-Scargle technique to high-cadence optical light curves from the Transiting Exoplanet Survey Satellite ( TESS ), we detected several significant oscillation modes in the 3 to 24 hour period range for each source. We also measured the IR excess (a proxy for disk growth) of those five sources, using J-band light curves from Palomar Gattini-IR. In four of the five sources, we found anti-correlations between the IR excess and the amplitude of the main flux oscillation modes. This result is inconsistent with the conventional idea that non-radial pulsations drive mass ejections. We propose an alternative scenario where internal temperature variations in the Be star cause transitions between pulsation-active and mass-ejection-active states.  \nKey words: stars: early-type - stars: emission-line, Be - stars: oscillations - stars:rotation - binaries: close - X-rays:binaries', '1 INTRODUCTION': 'A Be X-ray binary (BeXB) is a binary system composed of a compact object (a neutron star in most cases) and a Be star, an early type star with a circumstellar disk (Reig 2011). In addition to standard observational features of Be stars such as double-peaked emission lines and infrared (IR) excess (the latter usually explained as thermal bremsstrahlung from the disk), BeXBs exhibit various peculiar phenomena caused by the interaction of the compact object with the circumstellar disk. The most notable of them are normal (Type-I) and giant (Type-II) outbursts. Normal outbursts occur regularly and (quasi-) periodically, and are caused by the compact object capturing mass from the circumstellar disk when passing around the periastron. Their typical duration and luminosity are about 0.2-0.3 times the or-  \nbital period and ≲ 10 37 erg s -1 , respectively. On the other hand, giant outbursts are rarer, more luminous and unpredictable; their mechanism is not well understood. They last several times the orbital period, and their luminosity often exceeds ∼ 10 38 erg s -1 (Reig & Nespoli 2013), the Eddington limit of a neutron star. In at least three well-studied cases (SMC X-3: Tsygankov et al. 2017, Townsend et al. 2017; Swift J0243.6+6124: Wilson-Hodge et al. 2018, AlfonsoGarzón et al. 2024; RX J0209.6-7427: Vasilopoulos et al. 2020, Hou et al. 2022, West et al. 2024), the peak X-ray luminosity reached ≈ 2 × 10 39 erg s -1 , placing those BeXBs into the ultra-luminous Xray (ULX) pulsar class. Therefore, giant-outbursting BeXBs might be a significant component of the observed ULX population in the local universe (Earnshaw et al. 2018; Kuranov et al. 2020; Gúrpide et al. 2021; Misra et al. 2024). In particular, luminous BeXBs can be the dominant component of the X-ray binary population in galaxies where star formation has peaked ≈ 25-60 Myr ago (An-  \nl. 2010; Antoniou & Zezas 2016), as opposed to more recent (dominated by high-mass X-ray binaries) or older (dominated by low-mass X-ray binaries) star-forming environments. Thus, understanding the connection between donor-star variability and giant outbursts in BeXBs is an important step for the modelling of X-ray binary populations in external galaxies. Such studies also provide observational constraints to theoretical model of super-critical accretion onto magnetised neutron stars (Tsygankov et al. 2017; Mushtukov et al. 2017, 2019).  \nSeveral studies (Martin et al. 2011, 2014; Okazaki et al. 2013; Reig & Blinov 2018) suggest that giant outbursts happen when the outer radius of the circumstellar disk (possibly warped, inclined or eccentric) reaches the orbit of the compact object and the latter captures a large amount of mass in a short time. A circumstellar disk is formed via outwards viscous diffusion of mass ejected from the stellar surface (Rivinius et al. 2013; Rímulo et al. 2018). However, the mass ejection mechanism is still unclear. Although Be stars are generally rapid rotators, their rotation velocities are typically only about 70 per cent of their critical velocity (Rivinius et al. 2013). This means that their rotation alone cannot explain mass ejections. Nonradial pulsations are an alternative mechanism proposed to explain mass ejections (Bagnulo et al. 2012). In support of this scenario, observations of the B2Ve star 𝜇 Cen showed a coincidence of mass ejections with constructive interference of the pulsational velocity fields (Rivinius et al. 1998). However, theoretically, non-radial pulsations limited by the speed of sound seem to be insufficient for disk formation (Smith et al. 2012; Torrejón et al. 2012). Addressing this mass ejection problem is the main objective of this work.  \nThe operational launch of the Transiting Exoplanet Survey Satellite ( TESS ; Ricker et al. 2015) in 2018 was a notable event for the asteroseismology of massive stars, including Be stars. Although the main purpose of TESS is to explore exoplanets, its high-precision and high-cadence continuous light curves are very useful for the analysis of periodic flux oscillations of massive stars due to their pulsations and rotation. Several TESS studies of Be stars and BeXBs have already shed new light in this field. Balona & Ozuyar (2020) examined TESS light curves of 57 Be stars, and found that 74 per cent of their targets show a single or harmonic double peak in the 0.5-3 d -1 frequency range of their power spectra. They argued that such peaks are caused by stellar rotation. Ramsay et al. (2022) investigated TESS light curves of 23 high-mass X-ray binaries, and confirmed the detection of 0.1-1 d -1 oscillations in the power spectra of all their targets. Reig & Fabregat (2022) used TESS to study the power spectrum morphology of 22 BeXBs, and found that BeXBs and classical Be stars are indistinguishable in terms of pulsational characteristics.  \nIn this work, we analyse the TESS light curves of BeXBs to investigate long-term variations in the activity level of Be donor stars. We compare the results with the long-term multi-wavelength light curves in order to constrain the mass ejection mechanism of the Be donor star. This article is composed as follows. In Section 2, we provide a brief description of the target selection. Section 3 describes the observational data used in this work: TESS , Swift , MAXI, ZTF, ATLAS, and Gattini-IR. We explain our Lomb-Scargle analysis of the TESS light curves and the IR excess evaluation in Section 4. Section 5 shows long-term multi-wavelength light curves and amplitude spectra of the target sources. We highlight the discovery of an anticorrelation between the IR excess and the amplitude of the main flux oscillation modes. In Section 6, we discuss possible interpretations of this observed anti-correlation. We propose that internal temperature variations in the Be star cause transitions between pulsation-active and mass-ejection-active states. Finally, Section 7 concludes this article.', '2 TARGET SELECTION': 'For BeXBs which have an optical magnitude of ≲ 14 mag and available TESS light curves, and are in the northern sky, we checked multi-wavelength light curves (cf. Section 3.2) in the modified Julian date (MJD) range of 58400-60000 and found that five of them showed significant variations in X-ray, optical and near-infrared (NIR) band (cf. Section 5.1). In most cases, optical and infrared (OIR) variations of Be stars are derived by the growth or the decay of their circumstellar disks. Thus we selected these five BeXBs as the targets for this study. Table 1 describes the selected five sources: 4U0115+634, GRO J2058+42, RX J0440.9+4421, SAX J2103.5+4545, and V0332+53.  \nThese are already known to exhibit long-term multi-wavelength activities such as giant outbursts, disk-derived OIR variations. For 4U0115+634, long-term variations have been studied in Negueruela & Okazaki (2001) and Reig et al. (2007), and they reported quasiperiodic disk variations with an approximately 5-year cycle and giant outbursts linked to them. Furthermore, variations in polarization were observed that appear to be related to disturbances in the disk during the giant outburst (Reig & Blinov 2018). For GRO J2058+42, Reig et al. (2023) confirmed sinusoidal variations in optical light curves and the H 𝛼 equivalent width with a period of about 9.5 years and occurrences of giant outbursts related to them, and stated that the symmetry of the double-peaked emission lines suggested that the disk was not warped. Reig et al. (2005) reported long-term variability in OIR light curves and H 𝛼 equivalent width of RX J0440.9+4431, and suggested that disturbances in the disk due to interactions with neutron stars can only occur if the disk is sufficiently developed, based on V/R 1 variations. Long-term activities of SAX J2103.5+4545 are studied in Reig et al. (2010) and Camero et al. (2014). This source shows highly correlated X-ray and OIR variability, and X-ray flares linked to peaks of OIR light curves. One notable feature is that the orbital period of this source ( 𝑃 orb = 12 . 67 d; Camero Arranz et al. 2007) is particularly short for known BeXBs. Caballero-García et al. (2016) carried out a long-term X-ray and OIR observation of V0332+53, and suggested that the inner regions surrounding the magnetosphere was visualized during the lowest flux state based on the presence of X-ray quasi-periodic oscillation at that time.', '3 OBSERVATIONS': 'Weutilised high-cadence optical light curves of TESS , and long-term multi-wavelength light curves.', '3.1 Transiting Exoplanet Survey Satellite': 'Transiting Exoplanet Survey Satellite ( TESS ; Ricker et al. 2015) is a optical space telescope operated by National Aeronautics and Space Administration (NASA) and Massachusetts Institute of Technology (MIT), and its main purpose is to explore exoplanets by observing their transits. Full frame images (FFIs) at 10 or 30-minutes intervals are available from the Mikulski Archive for Space Telescopes (MAST).  \nWe carried out the photometry of FFIs with the pipeline of TESS Transients project (Fausnaugh et al. 2021, 2023). This pipeline evaluates a brightness of the star with difference imaging and fitting the point-spread function, and is expected to be more accurate than a simple aperture photometry for fainter sources ( TESS magnitude ≳ 12). FFIs record photo-electron count rate, and we converted it  \nTable 1. Target BeXBs.  \n- [a] Wenger et al. (2000); [b] Paegert et al. (2022); [c] Reig & Fabregat (2015); [1a] Tamura et al. (1992); [1b] Negueruela & Okazaki (2001);\n- [2a] Wilson et al. (2005); [2b] Kızıloˇglu et al. (2007); [3a] Ferrigno et al. (2013); [3b] Reig et al. (2005); [4a] Camero Arranz et al. (2007);\n- [4b] Reig et al. (2004); [5a] Negueruela et al. (1999); ∗ Distances are calculated from Gaia DR3 parallaxes (Gaia Collaboration et al. 2023).  \nFigure 1. Energy bands of multi-wavelength light curves. Data of ZTF and ATLAS filters were obtained from the SVO Filter Profile Service ( http://svo2.cab.inta-csic.es/svo/theory/fps/index.php ). As for Gattini-IR, it is drawn with the value of the Cousins J-band filter.  \n<!-- image -->  \nto energy flux density using the following relationship based on a documentation in an official web-page of TESS Transients project 2 :  \n𝐹 𝐼 = 1 . 61 × 10 -5 GLYPH<20> Jy counts s -1 GLYPH<21> , (1)  \nwhere 𝐼 is the photo-electron count rate, and 𝐹 is an energy flux density in the TESS pass-band. The interval of FFIs was changed from 30 to 10 minutes at the end of the prime mission in July 2020 (Ricker 2021). We unified the sampling rate of the light curves by binning at 30-minutes widths during the analysis process (cf. Section 4.3.1). In addition, we obtained TESS magnitude (Tmag) of target sources from TESS Input Catalog (Stassun et al. 2018, 2019).', '3.2 Long-term multi-wavelength light curves': 'Weutilised X-ray light curves of Swift /BAT and MAXI/GSC, optical light curves of ZTF and ATLAS, and NIR light curves of GattiniIR, covering MJD range of ∼ 58400-60000. Figure 1 shows energy bands of multi-wavelength light curves.', '3.2.1 Neil Gehrels Swift Observatory': 'Neil Gehrels Swift Observatory ( Swift ; Gehrels et al. 2004) is a gamma-ray burst survey satellite operated by NASA. We downloaded  \nthe Burst Alert Telescope (BAT) 15-50 keV daily light curves of target BeXBs, from Swift /BAT Hard X-ray Transient Monitor webpage (Krimm et al. 2013). Seven days binning was performed to improve signal-to-noise ratio for each light curves.', '3.2.2 Monitor of All-sky X-ray Image': '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).', '3.2.3 Zwicky Transient Facility': "Zwicky Transient Facility (ZTF; Bellm et al. 2019) is an optical widefield 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 𝑔 ' and 𝑟 ' -band light curves of the targets acquired from ZTF Public Data Release 15 (DR15; Masci et al. 2019), which covers the epoch of ∼ 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 ( ∼ 12.5 AB).", '3.2.4 Asteroid Terrestrial-impact Last Alert System': '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 𝑐 and 𝑜 -band light curves covering 58000-60000 MJD via ATLAS Forced Photometry web interface (Smith et al. 2020). The saturation magnitude of ALTAS is ∼ 12.5 AB, about same as ZTF. Thus, the light curve for RX J0440.9+4431, although it can be created, is not reliable.', '3.2.5 Gattini-IR': 'Gattini-IR (De et al. 2020) is a NIR wide-field survey project using 30 cm robotic 𝐽 -band telescope system at Palomar Observatory, operated by Caltech. Gattini-IR sweeps ∼ 7500 deg 2 of the sky at one night with an upper-limit magnitude of ∼ 15.7 (AB), and scans the northern-sky observable from Palomar (Dec. ≳ -30 · ) with a cadence of ∼ 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  \nTable 2. Coefficients of Equations (3) and (4) (Torres et al. 2010).  \ntransients. We obtained the 𝐽 -band light curves of the targets generated by GDPS. The data-set covers the epoch of ∼ 58400-59900 MJD.', '4 DATA ANALYSIS': 'We first estimated Be donor properties to evaluate the IR excess and the critical velocity. Next, the IR excess was evaluated to investigate the evolution of the disks. Finally, a flux periodicity was analysed applying the Lomb-Scargle technique (Lomb 1976; Scargle 1982) to TESS light curves. We utilised curve\\_fit function of scipy.optimize library (Virtanen et al. 2020) for fitting, and LombScargle class of astropy.timeseries library (Astropy Collaboration et al. 2013, 2018, 2022) as an implementation of the Lomb-Scargle technique.', '4.1 Estimation of Be donor properties': 'We estimated an effective temperature 𝑇 eff , radius 𝑅 ∗ , and mass 𝑀 ∗ for each Be donor of target BeXBs by fitting of an optical spectral energy distribution (SED) with a simple blackbody model. First, the optical SED of the target was created using catalogued optical magnitudes, a color excess 𝐸 ( 𝐵 -𝑉 ) , and a distance. The reason for not using the NIR data (e.g. 2MASS) is that it is more susceptible to the IR excess than the optical band and the blackbody radiation may not be dominant. We used the Cardelli et al. (1989) extinction function to derive extinctions for each band. This SED was then fitted with a model spectrum of a spherical blackbody to obtain 𝑇 eff and 𝑅 ∗ . Finally, 𝑀 ∗ was calculated using the empirical relations of Torres et al. (2010) with some modifications:  \n𝑋 = log GLYPH<18> 𝑇 eff 12 . 6 [ kK ] GLYPH<19> , (2)  \nlog GLYPH<18> 𝑀 ∗ 𝑀 ⊙ GLYPH<19> = 𝑎 1 + 𝑎 2 𝑋 + 𝑎 3 𝑋 2 + 𝑎 4 𝑋 3 + 𝑎 5 ( log 𝑔 ) 2 + 𝑎 6 ( log 𝑔 ) 3 ,  \n(3)  \nlog GLYPH<18> 𝑅 ∗ 𝑅 ⊙ GLYPH<19> = 𝑏 1 + 𝑏 2 𝑋 + 𝑏 3 𝑋 2 + 𝑏 4 𝑋 3 + 𝑏 5 ( log 𝑔 ) 2 + 𝑏 6 ( log 𝑔 ) 3 , (4)  \nwhere 𝑔 is a surface gravity normalised in unit of cm s -2 , 𝑀 ⊙ and 𝑅 ⊙ are the solar mass and radius, respectively, and 𝑎 1-6 and 𝑏 1-6 are coefficients. The base of the logarithm here is 10. Values of 𝑎 1-6 and 𝑏 1-6 are listed in Table 2. Specifically, we solved Equation (4) with a numerical technique to obtain log 𝑔 from 𝑇 eff and 𝑅 ∗ , and calculated 𝑀 ∗ by substituting 𝑇 eff and log 𝑔 into Equation (3). Note that we ignored metallicity [Fe/H] terms, which are in original equations, since all targets in this study are galactic sources. We obtained the optical multi-band magnitudes from APASS DR10 (Henden et al. 2018), PS1 DR1 (Flewelling et al. 2020), or Table 2 of Reig &  \nFabregat (2015). The indeterminacy of the result was estimated by the bootstrap re-sampling. If multiple results were obtained by multiple catalogues, we selected the one for which 𝑇 eff is most plausible to the typical value of its spectral type (cf. Table A1).  \nWeestimated a blackbody 𝐽 -band magnitude 𝑚 bb 𝐽 , critical velocity v b , and critical frequency 𝜈 b . 𝑚 bb 𝐽 was obtained by extrapolating the obtained blackbody SED model to the 𝐽 -band. This is the expected 𝐽 -band magnitude of the Be donor itself excluding the IR excess. 𝜈 b is a rotation frequency when spinning at v b . We calculated v b and 𝜈 b assuming a sphere:  \nv b = √︂ 𝐺𝑀 ∗ 𝑅 ∗ , (5) 𝜈 b = v b 𝜋𝑅 , (6)  \n2 ∗  \nwhere 𝐺 is the gravitational constant. Since the rotation of the Be star does not reach v b in most cases (Rivinius et al. 2013), 𝜈 b can be regarded as the upper limit of the rotational frequency. If the contribution of the IR excess to the catalog optical magnitudes is not negligible, it causes the color to be redder and the magnitude to be larger, and results in underestimation of 𝑇 eff and overestimation of 𝑅 ∗ and 𝑚 bb 𝐽 . Estimated properties are summarised in Table 3. At temperature 𝑇 ∼ 30000 K, the Planck distribution peak is at ultraviolet-band ( 𝜈 ∼ 10 15 Hz), and 𝑇 change produces only a small difference in slope in optical band. Therefore, the temperature estimation in this method is highly sensitive to small differences in the slope of the SED, resulting in a large error of 𝑇 eff due to a slight uncertainty of 𝐸 ( 𝐵 -𝑉 ) . Note that v b is consistent with the literature value of v sin 𝑖 (i.e. v b ≳ v sin 𝑖 ) for all sources (cf. Table 1).', '4.2 IR excess evaluation': 'We evaluated the IR excess as a measure of disk development. The 𝐽 -band flux observed by Gattini-IR was averaged per sector and converted into magnitudes, defined as 𝑚 𝐽 . We used 𝑚 𝐽 -𝑚 bb 𝐽 to evaluate the amount of the IR excess (Figure 2). The smaller the value of 𝑚 𝐽 -𝑚 bb 𝐽 , the greater the IR excess. The Gattini-IR magnitudes were systematically dimmer than 𝑚 bb 𝐽 in 4U0115+634. This can be qualitatively explained by overestimation of 𝑚 bb 𝐽 due to some factors such as underestimation of 𝑇 eff (cf. Section 4.1). In any case, the IR excess evaluated in this way can be used to discuss long-term variations of one source, but not for cross-source comparisons.', '4.3 Analysis of flux periodicity': 'We made and analyzed amplitude spectra of TESS light curves for each sector of each targets to investigate the flux periodicity. This analysis was carried out with the following procedure:  \n- (i) Data reduction to raw TESS light curves,\n- (ii) Amplitude spectrum estimation using the Lomb-Scargle technique,\n- (iii) Decomposition of the spectrum into peaks and the red noise.', '4.3.1 Data reduction': "The data reduction process consists of three stages: detrending, 5sigma 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 ≳ several days. For convenience, we call this process as 'detrending'. The width of the detrending window was set  \nTable 3. Estimated stellar properties of target BeXBs.  \n- ∗ Sources of the optical magnitudes used to derive the properties. [R&F] means Reig & Fabregat (2015).  \nFigure 2. IR excess evaluation of 4U0115+634. Gray square and blue circular markers correspond to the raw Gattini-IR 𝐽 -band light curve and the sectorby-sector averaged light curve, respectively. The horizontal dashed line is 𝑚 bb 𝐽 which is the estimated IR excess excluded magnitude, and the blue shaded regions represent sectors where TESS observed this source. We used the difference between the blue points and 𝑚 bb 𝐽 line as a measure of the IR excess. The fact that the light curve is sometimes systematically below 𝑚 bb 𝐽 (e.g. 58800-59500 MJD) may be due to an overestimation of 𝑚 bb 𝐽 .  \n<!-- image -->  \nto 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 ( ∼ 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.", '4.3.2 Amplitude spectrum': '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):  \n𝐴 ( 𝜈 ) = √︂ 4 𝑁 𝑃 ( 𝜈 ) , (7)  \nwhere 𝜈 is a frequency, 𝐴 ( 𝜈 ) is the amplitude spectrum, 𝑁 is a number of points in the light curve, and 𝑃 ( 𝜈 ) 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  \nFigure 3. Data reduction to the TESS light curve of V0332+53, sector 59. The three panels correspond to the 3-day detrending, 5-sigma clipping, and 30-min binning processes, respectively.  \n<!-- image -->  \nthe 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 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.', '4.3.3 Decomposition into peaks and red noise': "To extract the peaks from the spectrum, we estimated a red-noise component of the spectrum. First, we made a local 25th percentile  \nFigure 4. Dynamic amplitude spectrum of V0332+53, sector 59. The top panel is a light curve, and the bottom is a dynamic spectrum. The shaded area is where the frequency is less than 1/3 d -1 or the number of extracted light curve points is less than 160.  \n<!-- image -->  \nFigure 5. Averaged amplitude spectrum of V0332+53, sector 59. The top panel shows a 2D histogram of dynamic amplitude spectrum (Figure 4) projected along the time axis, and the bottom panel shows an averaged amplitude spectrum. The averaged amplitude spectrum can be paraphrased as the centroid of the histogram for each frequency. 'Log density' is the logarithm of density (arbitrary units) and its base is 10. The left shaded area is where the frequency is less than 1/3 d -1 . See also Section 4.3.3 for the description of the red noise.  \n<!-- image -->  \n(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-  \nFigure 6. Red noise estimation of V0332+53, sector 59. The smoothed local 25th percentile is the roughly estimated continuum, and white noise level is the constant component of the red noise model. 'Photometric Error' means a spectrum estimated from the photometric error. 'Normalized Amplitude' is the ratio of amplitude to the flux converted from Tmag. The left shaded area is where the frequency is less than 1/3 d -1 .  \n<!-- image -->  \nike term and a constant term:  \n𝐴 RN ( 𝜈 ) = √︄ 𝑎 2 0 1 + ( 𝜈 / 𝜈 0 ) 𝛾 + 𝐶 2 , (8)  \nwhere 𝜈 is a frequency, 𝐴 RN is a red noise amplitude spectrum, 𝑎 0 is a scaling factor, 𝜈 0 is a characteristic frequency, 𝛾 is an exponent, and 𝐶 is a white noise level. When 𝛾 = 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.  \nThe 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.", '5.1 Long-term multi-wavelengths activities': 'Figures 7 to 11 are long-term multi-wavelength light curves of five sources. They exhibited significant X-ray outbursts in both 2-20 keV and 15-50 keV light curves, and OIR variability of ≳ 0.5 mag. The X-ray light curves of SAX J2103.5+4545 during outbursts show two components: a sharp peak with a duration of about 10 days and a relatively gentle bump of about 100 days. This feature was also seen in previous outbursts of this source (cf. Camero et al. 2014).  \nThe Crab pulsar has a luminosity of ∼ 10 37 erg s -1 (Kirsch et al.  \nFigure 7. Multi-wavelength light curves of 4U0115+634. The gray shaded regions correspond to TESS observed sectors. An X-ray outburst can be confirmedat ∼ 58900MJD,whichisconsistent with the Atel report (Nakajima et al. 2020).  \n<!-- image -->  \nFigure 10. A similar figure to Figure 7 for SAX J2103.5+4545. Recurrent X-ray/OIR outbursts can be confirmed at ∼ 58600, 59000, and 59500 MJD.  \n<!-- image -->  \nFigure 8. Asimilar figure to Figure 7 for GRO J2058+42. An X-ray outburst can be confirmed at ∼ 58500 MJD, which is consistent with the GCNC report (Barthelmy et al. 2019).  \n<!-- image -->  \n2005), a distance of ∼ 2 kpc (Lin et al. 2023), and is detected by the GSContheorder of 1 photons s -1 cm -2 (Morii et al. 2011). By comparing GSC count rates and distances of five targets with those of the Crab pulsar, the peak luminosities of these outbursts can be roughly estimated to be ≳ 10 37 erg s -1 . In addition, the orbital periodicity expected for normal outbursts was not confirmed, and the duration of outbursts was ∼ 100 days, which is longer than orbital periods other than RX J0440.9+4431. Outbursts of SAX J2103.5+4545 appear to be periodic, but their interval ( ∼ 200 days) is not consistent with its orbital period (12.67 days). Therefore, all detected X-ray outbursts can be classified as giant outbursts.  \nThe time scales for all OIR light curve variations were about hundreds of days, which are consistent with the typical time scale of the circumstellar disk evolution. Figure 12 shows the OIR color variations of 4U0115+634. The spectral index shown in the figure represents the color, and is calculated for given two bands as follows:  \n𝛼 ( 𝜈 0 , 𝜈 1 ) = log GLYPH<0> 𝑓 𝜈 0 / 𝑓 𝜈 1 GLYPH<1> log ( 𝜈 0 / 𝜈 1 ) , (9)  \nwhere 𝜈 0 and 𝜈 1 , and 𝑓 𝜈 0 and 𝑓 𝜈 1 are photon frequencies and fluxes in the given two bands, respectively, and 𝛼 is a spectral index. This corresponds to the slope of the line connecting two points in the dou-  \nFigure 9. A similar figure to Figure 7 for RX J0440.9+4431. Optical light curves exceeded the saturation level and are unavailable for analysis. An X-ray outburst can be confirmed at ∼ 60000 MJD, which is consistent with the Atel report (Nakajima et al. 2022).  \n<!-- image -->  \nFigure 11. Asimilar figure to Figure 7 for V0332+53. An X-ray outburst can be confirmed at ∼ 58700 MJD.  \n<!-- image -->  \nFigure 12. OIR color variations of 4U0115+634. Top and bottom panel show the J-band light curve and OIR color variations, respectively. The vertical axis in the bottom panel is the spectral index, which is redder at the bottom and bluer at the top.  \n<!-- image -->  \nand a larger value indicates a bluer spectrum. The figure shows a "redder when brighter" tendency that are common for OIR variations of Be stars. This tendency was confirmed for all targets except RX J0440.9+4431 whose optical light curves were unavailable. Therefore, all of the observed long-term OIR variations can be attributed to the disk evolution. Thus, observed 𝐽 -band flux, and hence 𝑚 𝐽 -𝑚 bb 𝐽 , can be used as a measure of the disk development.  \nRecent studies on giant outbursts (Martin et al. 2011, 2014; Okazaki et al. 2013) suggested that such events should occur when a neutron star passes near a sufficiently developed disk that extends to the orbit of the neutron star. Since the orbital period of BeXBs is typically shorter than the time scale of the disk growth, giant outbursts are expected to occur at or shortly before the maximum IR excess. However, all detected giant outbursts were initiated months to hundreds of days after the peak of the 𝐽 -band light curve. The same was observed in the 2017 giant outburst of Swift J0243.6+6124 (AlfonsoGarzón et al. 2024). This discrepancy can be explained by assuming that the peak IR excess and the maximum outer disk radius are not  \nTable 4. Detected peak frequencies of five sources. Peaks with a check mark in the \'Harm.\' column are those which have a possible integer multiple relationship to the frequency of the other peaks.  \nsimultaneous. Since the mass is supplied to the disk from the Be star, i.e., from the innermost edge of the disk, there should be a time lag of about the viscous time scale between the mass supply and the outward extension of the disk. In other words, the total mass and outermost radius never reach the maximum at the same time, and the former does earlier. In addition, since the volume emissivity of thermal bremsstrahlung is proportional to the square of the electron density, IR excess should depend on both disk radius and mass, and the peak of IR excess should also be before the radius peak. Hence, the observed delay can be attributed to the time lag between mass supply and disk extension.', '5.2 Amplitude spectra': 'We summarised detected peak frequencies in Table 4, and showed red noise subtracted amplitude spectra of five sources in Figures 13 to 17. When the ratio of two peak frequencies matched an integer within the 1 𝜎 error range, we determined that the pair of peaks has possible harmonic relationship, except in cases where the error of ratio exceeded 0.5. The possible harmonic relationships were confirmed for only two of the five, 4U0115+634 and RX J0440.9+4431. These were the only two having peaks of 𝜈 ≲ 𝜈 b . In all five sources, the maximum peak was detected at 𝜈 = 2-4 d -1 with a normalized amplitude of ∼ 1 per cent. There was at least one peak per source where the normalized amplitude varied clearly between sectors. And the maximum peak was the amplitude varying peak for all sources.', '5.3 Peak amplitude and the IR excess correlations': "Figures 18 to 22 show the relationship between peak amplitudes and the IR excess for five sources. For all sources except GRO J2058+42, at least one peak exhibited significant anti-correlation between its amplitude and the IR excess. Other peaks did not show significant correlations. We refer these peaks as 'anti-correlated peaks' and 'non-correlated peaks', respectively. Except for GRO J2058+42, the maximum peak belonged to the anti-correlated peaks. The ratio of the maximum to minimum amplitudes of the anti-correlated peaks  \nFigure 13. Red noise subtracted amplitude spectra of 4U0115+634. 'Normalized Amplitude' is the amplitude normalised by the flux calculated from the Tmag of the source. 'Red Noise × 2' line is the threshold of the peak detection. Vertical dashed lines indicate positions of detected peaks.  \n<!-- image -->  \nFigure 14. A similar figure to Figure 13 for GRO J2058+42.  \n<!-- image -->  \nwas typically 5-6. For RX J0440.9+4431, both the anti-correlated and non-correlated peaks coexist (e.g. peaks at 𝜈 = 2 . 2 , 2 . 6).", '6.1 Origin of periodic flux oscillations': 'The origin of the periodic flux oscillation, which appears as a peak in the amplitude spectrum, is the rotation or pulsation of the Be donor star. According to Balona (2019) and Balona & Ozuyar (2020), rotation-derived flux oscillations in massive stars are originated from gas clouds co-rotating with the star, and these appear in the power spectrum as clusters of peaks or broad humps, and integer multiplied harmonics. Furthermore, the frequency of the rotation-derived fundamental mode should be smaller than 𝜈 b , because the rotation of the Be star does not reach a critical speed in most cases. In light of this, peaks of 𝜈 < 2 d -1 in 4U0115+634 and RX J0440.9+4431 may be  \nFigure 15. A similar figure to Figure 13 for RX J0440.9+4431.  \n<!-- image -->  \nFigure 16. A similar figure to Figure 13 for SAX J2103.5+4545.  \n<!-- image -->  \nFigure 17. A similar figure to Figure 13 for V0332+53.  \n<!-- image -->  \nFigure 18. IR excess and peak amplitude relationship in 4U0115+634. (top) Red noise subtracted normalised amplitude spectrum. Basically the same as Figure 13. (bottom) The further to the right, the larger the IR excess. Points with the same marker style and connected by dashed lines correspond to the peak of the same frequency and are plotted with values in multiple sectors. The color of the marker is the same as that of the vertical line in the top panel. In the upper right box are the correlation coefficients for each peak. Note that dashed lines simply connect adjacent points and ignore the time series.  \n<!-- image -->  \nFigure 19. A similar figure to Figure 18 for GRO J2058+42.  \n<!-- image -->  \nFigure 20. A similar figure to Figure 18 for RXJ 0440.9+4431.  \n<!-- image -->  \nFigure 21. A similar figure to Figure 18 for SAX J2103.5+4545.  \n<!-- image -->  \nrotation-derived fundamental modes, and their integer multiples may be harmonics. On the other hand, the other peaks do not have rotationderived features, and are considered to be pulsation-derived. However, the possibility that they are harmonics of the rotation-derived peaks cannot be completely ruled out. Furthermore, it should be noted that the ratio of the frequencies of the rotation-derived and pulsation-derived peaks could be close to an integer by chance, or that either the fundamental or harmonics could not be detected due to lack of signal-to-noise ratio.  \nFigure 22. A similar figure to Figure 18 for V0332+53.  \n<!-- image -->', '6.2 Interpretation of observed anti-correlations': 'The anti-correlation between flux oscillation amplitude and the IR excess is a notable result. According to the conventional idea of the pulsation-driven mass ejection, the pulsations should be active when the disk is growing, and thus the flux oscillation amplitude should be positively correlated with the IR excess, but the result was the opposite. In other words, if the anti-correlated peak amplitude variation we found is derived by the variation in the pulsation amplitude, it may disprove pulsation-driven mass ejection. However, if the pulsation amplitude did not change but only the flux oscillation amplitude changed, or if the anti-correlated peak is not pulsation-derived but rotation-derived, it is possible to explain the result without denying the pulsation-driven mass ejection. We examined following four scenarios to interpret the result:  \n- 1st: Co-rotating gas cloud scenario,\n- 2nd: Fully covered photosphere scenario,\n- 3rd: Partially covered photosphere scenario,\n- 4th: Internal state transition scenario.  \nWe have to mention the possibility that peak amplitude variations were caused by low-frequency beats ( 𝜈 ∼ 1 / 100 d -1 ), rather than the astrophysical mechanism. However, in this case, the anti-correlation can only be explained as being due to observations made at phases that happen to be so, which is inconsistent with the fact that anticorrelations were confirmed in four of the five sources. It is also possible that the anti-correlations are outliers. Our data set has too small a sample size to examine statistical significance.', '6.2.1 1st: Co-rotating gas cloud scenario': '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  \nFigure 23. Schematics of the first scenario. These illustrations view the Be star along its rotation axis. On the left is the case of a low mass ejection, which has a larger dynamic range of co-rotating gas cloud density along the rotational direction. On the other hand, the right side is when mass ejection is active, and the density distribution in the rotational direction is homogenised by the increased amount of co-rotating gas, resulting in a smaller dynamic range of density.  \n<!-- image -->  \ndetermined 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 noncorrelated 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.', '6.2.2 2nd: Fully covered photosphere scenario': '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, weexpect to observe the peak amplitude decrease in conjunction with the disk growth. This logic can also be applied to explain rotationderived 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.  \nBecausethesituation in this scenario is simple, it is easy to estimate the density of the circumstellar materials. The following relationship  \nFigure 24. Schematics of the second scenario. The green arrows are photon streams from the Be star, and the size indicates the number of photons pass through without scattering. In this scenario, the photons coming from the Be star are uniformly scattered by surrounding circumstellar materials at the same rate.  \n<!-- image -->  \nis expected to be satisfied when photons pass through a material: 𝜏 = e -𝑁𝜎 , (10)  \nwhere 𝜏 is a transmittance, 𝑁 is a electron column density, and 𝜎 is a scattering cross section. Substituting 𝜏 = 0 . 2 and the Thomson scattering cross section 𝜎 𝑇 = 6 . 7 × 10 -25 cm 2 for 𝜎 , the column density is calculated as 𝑁 = 2 . 4 × 10 24 cm -2 . Assuming a spherical symmetric stellar wind as scattering materials, an electron number density 𝑛 is expected to follow an inverse square law with respect to a distance from the star 𝑟 due to the conservation of particle number, when ignoring particle acceleration:  \n𝑛 = 𝑛 0 GLYPH<18> 𝑅 ∗ 𝑟 GLYPH<19> 2 , (11)  \nwhere 𝑛 0 is a number density at the surface and 𝑅 ∗ is a radius of the star. 𝑁 is then equal to integral of 𝑛 from the stellar surface to infinity:  \n𝑁 = ∫ 𝑛 d 𝑟  \n∞ 𝑅 ∗ (12)  \n= 𝑛 0 𝑅 ∗ . (13)  \nIf 𝑅 ∗ = 10 𝑅 ⊙ , then 𝑛 0 = 3 . 4 × 10 12 cm -3 . Assuming that the material is hydrogen, the numbers of electrons and nucleons are equal, thus a density at the surface is 𝜌 = 6 × 10 -12 g cm -3 . This is comparable to the observed density of the inner disk ( = 10 -13 -10 -10 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 ∼ 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  \nFigure 25. Roughly estimated density profile. The unit of density is g cm -3 and the base of logarithm is 10. The flatness in 𝑢 ∼ -1 is due to the neglect of particle acceleration near the stellar surface. The range of the vertical axis is the same as in Figure 15 of Curé (2004).  \n<!-- image -->  \nas 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 ∼ 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.  \nIncidentally, 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.', '6.2.3 3rd: Partially covered photosphere scenario': '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 har-  \nFigure 26. Schematics of the third scenario. In this scenario, the flux oscillation amplitude is reduced because the region near the Be star or the photosphere contributing to the flux oscillation is covered by something around the Be star (in this figure, the disk) when viewed from our perspective. If the flux oscillations are caused by a larger region, this scenario cannot account for the amplitude variation of ∼ 80 percent.  \n<!-- image -->  \nmonics. The local brightness at the expanding region is expected to be smaller due to gravitational darkening. In low𝑙 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𝑙 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𝑙 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.', '6.2.4 4th: Internal state transition scenario': '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 𝜅 -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 𝜅 is positively correlated with the temperature 𝑇 , 𝜅 variations provide positive feedback to 𝑇 variations, causing growth of thermodynamic fluctuations and thus pulsations. Since the amount of this feedback should depend on the  \nFigure 27. Schematic 𝑇 -𝜅 curve and temperature dependence of d 𝜅 / d 𝑇 . Pulsations should be most active where d 𝜅 / d 𝑇 is largest. And it is qualitatively expected that the pulsation will be inactive at higher or lower temperature regions.  \n<!-- image -->  \nslope of the 𝑇 -𝜅 curve (i.e. d 𝜅 / d 𝑇 ), the intensity of 𝜅 -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 𝜅 -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 𝑇 -𝜅 curve. And therefore, d 𝜅 / d 𝑇 becomes smaller in the higher temperature as it is closer to the Z-bump peak. Conversely, moving to the low temperature side also reduces d 𝜅 / d 𝑇 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.  \nSince 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.', '7 CONCLUSIONS': "The BeXBs exhibit characteristic and poorly understood phenomena such as two types of outbursts. To clarify them, we need to understand the nature of the circumstellar disk of the Be donor star. In particular,  \nTable 5. Summary of four scenarios. It shows whether each scenario can explain the observed results. ' ⃝ ', ' × ', and ' △ ' indicate explainable, inexplicable, and conditionally explainable, respectively.  \nthe mechanism driving the mass ejection of Be stars, which grows their circumstellar disks, is one of the critical unresolved issues.  \nWe analysed long-term variations in the flux periodicity of five galactic BeXBs using TESS light curves, and compared them with the long-term multi-wavelength light curves, to constrain the behaviour of their circumstellar disk. We obtained X-ray, optical, and NIR light curves for Swift , MAXI, ZTF, ATLAS, and Gattini-IR. Then we confirmed that there was OIR variability on time scale of hundreds of days derived from the circumstellar disk and they exhibited giant outbursts. We extrapolated the 𝐽 -band magnitudes from the optical catalogued magnitudes using a simple blackbody model and subtracted them from the Gattini-IR 𝐽 -band light curve to evaluate the IR excess. In addition, we made amplitude spectra from TESS light curves, and extract peaks of periodic flux oscillations derived by the rotation or pulsations.  \nAs a result, we found the anti-correlations between peak amplitudes and the IR excess. We examined the four scenarios to explain these results. Table 5 shows whether each scenario can explain the observed results. The first scenario, which assumes that all pulsationderived oscillations appears as non-correlated peaks, is unreasonable, and the second and third scenarios, which attempt to explain anticorrelations by reprocessing of the circumstellar materials, have severe conditions to explain the results. Therefore, we determined that it is most reasonable to interpret the observed anti-correlations as indicating an true anti-correlation between pulsation amplitudes and disk growth. In other words, the conventional idea of pulsation-driven mass ejection of Be stars may be incorrect. From this perspective, the most plausible scenario is the fourth one that assumed Be star internal state transition.  \nThere were only two or three sectors with both TESS and Gattini-IR data available for each source, so it is difficult to examine the statistical significance of the confirmed anti-correlations. In addition, there were no TESS observations covering the near-infrared brightening phase when the mass ejection should be most active, for the target sources of this study. Therefore, further accumulation of data from future TESS observations is desirable. It is also important to check whether other BeXBs exhibit the same anti-correlation.", 'ACKNOWLEDGEMENTS': "M.Niwanowassupported by Grant-in-Aid for JSPS Fellows. This research was supported by JSPS KAKENHI Grant Number 23KJ0913, and 23K25910. This work includes data collected by the TESS mission, and funding for the TESS mission is provided by NASA's  \nScience Mission Directorate. This research has made use of the MAXI data provided by RIKEN, JAXA, and the MAXI team. We acknowledge the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project, supported by the National Science Foundation under Grant No. AST1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron, and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen's University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile. This paper makes use of data from the AAVSO Photometric All Sky Survey, funded by the Robert Martin Ayers Sciences Fund and NSF (AST-1412587). The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia ( https://www.cosmos.esa.int/gaia ), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https: //www.cosmos.esa.int/web/gaia/dpac/consortium ). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work utilized Python libraries for analysis and visualization: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), pandas (Wes McKinney 2010), Matplotlib (Hunter 2007). This work made use of Astropy: 3 a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2013, 2018, 2022).", 'DATA AVAILABILITY': 'Swift /BATandZTFlightcurves used in this work are available in web pages of Swift /BAT Hard X-ray Transient Monitor 4 and NASA/IPAC Infrared Science Archive 5 , respectively. ZTF data we used are included in the ZTF public data release 15. The MAXI and ATLAS light curves in this work were generated by publicly available web services: MAXI on-demand 6 and ATLAS Forced Photometry 7 , respectively. As for the TESS light curves, the FFIs are available from MAST 8 . The Gattini-IR light curve is not publicly available.', 'REFERENCES': 'Stassun K. G., et al., 2018, AJ, 156, 102  \nTable A1. Typical properties of O-B type main sequence stars (Pecaut & Mamajek 2013). Extracted from https://www.pas.rochester.edu/ ~emamajek/EEM\\_dwarf\\_UBVIJHK\\_colors\\_Teff.txt .  \nStassun K. G., et al., 2019, AJ, 158, 138  \nTamura K., Tsunemi H., Kitamoto S., Hayashida K., Nagase F., 1992, ApJ, 389, 676  \nTonry J. L., et al., 2018, PASP, 130, 064505  \nTorrejón J. M., Schulz N. S., Nowak M. A., 2012, ApJ, 750, 75  \nTorres G., Andersen J., Giménez A., 2010, A&ARv, 18, 67  \nTownsend L. J., Kennea J. A., Coe M. J., McBride V. A., Buckley D. A. H., Evans P. A., Udalski A., 2017, MNRAS, 471, 3878  \nTsygankov S. S., Doroshenko V., Lutovinov A. A., Mushtukov A. A., Poutanen  \nJ., 2017, A&A, 605, A39  \nVasilopoulos G., et al., 2020, MNRAS, 494, 5350  \nVirtanen P., et al., 2020, Nature Methods, 17, 261  \nWenger M., et al., 2000, A&AS, 143, 9  \nWes McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Pro- ceedings of the 9th Python in Science Conference. pp 56 - 61, doi:10.25080/Majora-92bf1922-00a  \nWest B. F., Becker P. A., Vasilopoulos G., 2024, ApJ, 966, L5  \nWilson-Hodge C. A., et al., 2018, ApJ, 863, 9  \nWilson C. A., Weisskopf M. C., Finger M. H., Coe M. J., Greiner J., Reig P., Papamastorakis G., 2005, ApJ, 622, 1024  \nŠtefl S., Le Bouquin J. B., Carciofi A. C., Rivinius T., Baade D., Rantakyrö  \nF., 2012, A&A, 540, A76', 'APPENDIX A: TYPICAL PROPERTIES OF O-B TYPE MAIN-SEQUENCE STARS': 'Table A1 summaries typical properties of O-B type main sequence stars (Pecaut & Mamajek 2013).', 'APPENDIX B: A FRACTION OF PHOTONS RETURNING TO THE BE STAR DUE TO SCATTERING IN SPHERICALLY SYMMETRIC STELLAR WIND': 'A variation in the number of photons N as they pass through matter over a microscopic distance d 𝑟 can be written as follows:  \nd N = -𝑛𝜎 N d 𝑟, (B1)  \nwhere 𝑛 is a electron number density and 𝜎 is a scattering cross section. By solving Equation (B1) with substituting Equation (11), a  \nnumber of pass-through photons N pass is obtained:  \nN pass = N 0 exp GLYPH<20> -𝑛 0 𝜎𝑅 ∗ GLYPH<18> 1 -𝑅 ∗ 𝑟 GLYPH<19> GLYPH<21> , (B2) = N 0 exp h -𝑁𝜎 GLYPH<16> 1 -˜ 𝑟 -1 GLYPH<17> i , ( 𝑁 = 𝑛 0 𝑅 ∗ , ˜ 𝑟 = 𝑟 / 𝑅 ∗ ) (B3) = N 0 exp h -𝐴 GLYPH<16> 1 -˜ 𝑟 -1 GLYPH<17> i , ( 𝐴 ≡ 𝑁𝜎 )  \n(B4)  \nwhere 𝑅 ∗ is a stellar radius, N 0 and 𝑛 0 are a number of photons and an electron number density at the stellar surface ( 𝑟 = 𝑅 ∗ ), 𝑁 is an electron column density, and ˜ 𝑟 is a normalised radius. Note that Equation (B3) is consistent with Equation (10) in the limit of ˜ 𝑟 → ∞ . In terms of the model in Section 6.2.2, N pass corresponds to the number of photons passing a sphere of radius ˜ 𝑟 that have never experienced scattering. Since scattering experienced photons also pass through the same plane, N is expected to be between N pass and N 0 . The differential number of scattered photons is equal to the absolute value of the right side of Equation (B1):  \nd N scatt = 𝑛𝜎 N d 𝑟, (B5)  \n= 𝐴 N d˜ 𝑟 ˜ 𝑟 2 , (B6)  \nwhered N scatt is a number of photons scattered in a thin spherical shell of radius ˜ 𝑟 and thickness d˜ 𝑟 . If the scattering process is Thomson scattering and incident photons are unpolarized, the scattering is isotropic. In this case, the fraction of photons returning to the Be star is proportional to the solid angle of the star viewed from the scattering point. Assuming the spherical star, the solid angle of the star Ω can be calculated as follows:  \nΩ = ∫ sin -1 1 ˜ 𝑟 2 𝜋 sin 𝜃 d 𝜃,  \n0 (B7) = 2 𝜋 GLYPH<16> 1 -√︁ 1 -˜ 𝑟 -2 GLYPH<17> , (B8) = 2 𝜋𝐹 ( ˜ 𝑟 -1 ) . GLYPH<16> 𝐹 ( 𝑥 ) ≡ 1 -√︁ 1 -𝑥 2 GLYPH<17> (B9)  \nThen, the differential number of returning photons d N ret can be calculated from d N scatt and Ω :  \nd N ret = Ω 4 𝜋 d N scatt , (B10) = 𝐴 2 N 𝐹 ( ˜ 𝑟 -1 ) d˜ 𝑟 ˜ 𝑟 2 . (B11)  \nWe can obtain the total number of returning photons N ret by integrating d N ret from the stellar surface to the infinity:  \nN ret = ∫ ∞ 1 𝐴 2 N 𝐹 ( ˜ 𝑟 -1 ) d˜ 𝑟 ˜ 𝑟 2 . (B12)  \nIf N = N 0 , the fraction of returning photons N ret /N 0 is calculated as follows:  \nN ret N 0 = 𝐴 2 ∫ ∞ 1 𝐹 ( ˜ 𝑟 -1 ) d 𝑟 ˜ 𝑟 2 , (B13) = 𝐴 2 ∫ 1 0 𝐹 ( 𝑢 ) d 𝑢, GLYPH<16> 𝑢 = ˜ 𝑟 -1 GLYPH<17> (B14) = 𝐴 2 GLYPH<16> 1 -𝜋 4 GLYPH<17> = 0 . 11 𝐴. (B15)  \nSimilarly, the case N = N pass is as follows:  \nN ret N 0 = 𝐴 2 ∫ ∞ 1 exp h -𝐴 GLYPH<16> 1 -˜ 𝑟 -1 GLYPH<17> i 𝐹 ( ˜ 𝑟 -1 ) d˜ 𝑟 ˜ 𝑟 2 , (B16)  \n= 𝐴 2 e -𝐴 ∫ 1 0 𝐹 ( 𝑢 ) e 𝐴𝑢 d 𝑢, (B17)  \n= 𝐴 2 e -𝐴 GLYPH<18> e 𝐴 -1 𝐴 -∫ 1 0 √︁ 1 -𝑢 2 e 𝐴𝑢 d 𝑢 GLYPH<19> . (B18)  \nHere the first order modified Bessel and Struve function of the first kind, 𝐼 1 and 𝐿 1 can be written in following integral representations:  \n𝐼 1 ( 𝑥 ) = 2 𝜋 𝑥 ∫ 1 0 √︁ 1 -𝑢 2 cosh ( 𝑥𝑢 ) d 𝑢, (B19)  \n𝐿 1 ( 𝑥 ) = 2 𝜋 𝑥 ∫ 1 0 √︁ 1 -𝑢 2 sinh ( 𝑥𝑢 ) d 𝑢. (B20)  \nHence,  \n𝜋 𝐼 1 ( 𝑥 ) + 𝐿 1 ( 𝑥 ) 2 𝑥 = ∫ 1 0 √︁ 1 -𝑢 2 [ cosh ( 𝑥𝑢 ) + sinh ( 𝑥𝑢 )] d 𝑢, (B21) = ∫ 1 0 √︁ 1 -𝑢 2 e 𝑥𝑢 d 𝑢. (B22)  \nTherefore,  \nRight side of (B18) = 𝐴 2 e -𝐴 GLYPH<20> e 𝐴 -1 𝐴 -𝜋 𝐼 1 ( 𝐴 ) + 𝐿 1 ( 𝐴 ) 2 𝐴 GLYPH<21> , (B23)  \n= 1 2 GLYPH<26> 1 -e -𝐴 GLYPH<20> 1 + 𝜋 𝐼 1 ( 𝐴 ) + 𝐿 1 ( 𝐴 ) 2 GLYPH<21> GLYPH<27> . (B24)  \nSubstituting 𝜎 = 𝜎 T = 6 . 7 × 10 -25 cm 2 and 𝑁 = 2 . 4 × 10 24 cm -2 to Equations (B15) and (B24), the fraction is calculated as N ret /N 0 = 0 . 13 , 0 . 17. In other words, about 10-20 % of the total photons return to the Be star in the model considered in Section 6.2.2.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}