Gammaray binary LS I61 303 consists of a neutron star orbiting around a Be star with a period of P SUBorbSUB  26.5 days. Apart from orbital modulations the binary shows longterm flux variations with a superorbital period of inlineformula mmlmath overflowscrollmmlmsubmmlmrowmmlmiPmmlmimmlmrowmmlmrowmmlmisupmmlmimmlmrowmmlmsubmmlmommlmommlmn4.6mmlmnmmlmspace width0.25emmmlmspacemmlmrowmmlmrowmmlmi mathvariantnormalymmlmimmlmrowmmlmrowmmlmi mathvariantnormalrmmlmimmlmrowmmlmrowmmlmath inlineformula as seen in nearly all wavelengths. The origin of this superorbital modulation is still not well understood. Under the pulsar windstellar outflow interaction scenario we propose that the superorbital modulations of LS I61 303 could be caused by the precession of the Be disk. Assuming Xrays arise from synchrotron radiation of the intrabinary shock we develop an analytical model to calculate expected flux modulations over the orbital and superorbital phases. The asymmetric twopeak profiles in orbital light curves and sinusoidallike longterm modulations are reproduced under the precessing disk scenario. The observed orbital phase drifting of the Xray peak and our fitting of longterm Xray data indicate that the neutron star is likely orbiting around the star with a small eccentricity and periastron phase around SUBpSUB  0.6. We compare the Corbet diagrams of LS I61 303 with other BeXray binaries and the linear correlation in the inlineformula mmlmath overflowscrollmmlmsubmmlmrowmmlmiPmmlmimmlmrowmmlmrowmmlmisupmmlmimmlmrowmmlmsubmmlmommlmommlmsubmmlmrowmmlmiPmmlmimmlmrowmmlmrowmmlmiorbmmlmimmlmrowmmlmsubmmlmath inlineformula diagram suggests that the precession of the Be disk in LS I61 303 is induced by the tidal torque of its neutron star companion.
['Gamma-ray sources', 'X-ray binary stars', 'Compact binary stars', 'X-ray sources', '633', '1811', '283', '1822', 'Astrophysics - High Energy Astrophysical Phenomena']
A Precessing Stellar Disk Model for Superorbital Modulations of the GammaRay Binary LS I61 303
{'A precessing stellar disk model for superorbital modulations of the gamma-ray binary LS I+61 · 303': 'A. M. Chen, 1, 2 J. Takata, 3 and Y. W. Yu 4, 2  \n1  \nTsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 201210, China. Key Laboratory of Quark and Lepton Physics (MOE), Central China Normal University, Wuhan 430079, China.  \n- 2  \n3 Department of Astronomy, School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China. 4 Institute of Astrophysics, Central China Normal University, Wuhan 430079, China.', 'ABSTRACT': 'Gamma-ray binary LS I+61 · 303 consists of a neutron star orbiting around a Be star with a period of P orb /similarequal 26 . 5 d. Apart from orbital modulations, the binary shows long-term flux variations with a superorbital period of P sup /similarequal 4 . 6 yrs as seen in nearly all wavelengths. The origin of this superorbital modulation is still not well understood. Under the pulsar wind-stellar outflow interaction scenario, we propose that the superorbital modulations of LS I+61 · 303 could be caused by the precession of the Be disk. Assuming X-rays arise from synchrotron radiation of the intrabinary shock, we develop an analytical model to calculate expected flux modulations over the orbital and superorbital phases. The asymmetric two-peak profiles in orbital light curves and sinusoidal-like long-term modulations are reproduced under the precessing disk scenario. The observed orbital phase drifting of X-ray peak and our fitting of long-term X-ray data indicate that the neutron star is likely orbiting around the star with a small eccentricity and periastron phase around Φ p ∼ 0 . 6. We compare the Corbet diagrams of LS I+61 · 303 with other Be/X-ray binaries, and the linear correlation in the P sup -P orb diagram suggests that the precession of the Be disk in LS I+61 · 303 is induced by the tidal torque of its neutron star companion.  \nKeywords: Compact binary stars (283) - Gamma-ray sources (633) - X-ray sources (1822) - X-ray binary stars (1811)', '1. INTRODUCTION': "LS I+61 · 303 is a massive γ -ray binary that contains a compact object (CO) and a rapidly rotating B0Ve star (see Marcote 2017 for a review of this binary). Due to the unclear nature of CO, LS I+61 · 303 had been widely believed to be a microquasar for a long time (e.g.Gregory & Neish 2002; Bosch-Ramon & Paredes 2004; Gupta & Bottcher 2006; Massi, Migliari, & Chernyakova 2017; Jaron 2021), although some studies suggest that the CO is an energetic pulsar (Maraschi & Treves 1981; Lipunov & Nazin 1994; Dubus 2006; Dhawan, Mioduszewski, & Rupen 2006; Zdziarski, Neronov, & Chernyakova 2010). Recently, Weng et al. (2022) reported detections of transient radio pulses from LS I+61 · 303 with a period of 269 ms, which strongly suggests that the CO is a rotation-powered neutron star (NS) 1 . It also implies that the broadband emissions from the binary are likely resulting from the interaction between the relativistic pulsar wind with stellar outflows, similar to other γ -ray binaries like PSR B1259-63/LS 2883 and PSR J2032+4127/MT 91 213 (Tavani & Arons 1997; Takata et al. 2012, 2017; Petropoulou et al. 2018; Chen et al. 2019; Chen & Takata 2022).  \nSince its discovery by Gregory & Taylor (1978), LS I+61 · 303 has been intensively monitored by many telescopes across the whole electromagnetic spectrum (e.g., Hutchings & Crampton 1981; Taylor & Gregory 1982, 1984; Paredes et al. 1997, 2007; Albert et al. 2006, 2008; Chernyakova, Neronov, & Walter 2006; Chernyakova et al. 2012, 2023; Acciari et al. 2008, 2009, 2011; Hadasch et al. 2012; Aleksi'c et al. 2012; Jaron & Massi 2014; Archambault et al. 2016; Marcote et al. 2016; Zhang 2024). The system shows orbital modulated radiation from radio to very-high-energy (VHE) γ -rays with a period of P orb = 26 . 496 d. Anderhub et al. (2009) conducted a simultaneous multiwavelength  \nobservation of the binary, and they found a significant flux correlation between X-rays and VHE γ -rays. Orbital light curves in these two bands show a prominent peak around phase 0.6-0.7 and a second broader peak around phase 0.8-1.1. Chen & Takata (2022) proposed that this double-peak profile is due to interactions of the pulsar wind with the inclined Be disk.  \nIn addition to orbital modulations, LS I+61 · 303 also shows superorbital variations with a period of P sup /similarequal 4 . 6 yr as seen at nearly all wavelengths, including radio (Gregory 2002; Marcote et al. 2016; Jaron, Kiehlmann, & Readhead 2024), optical (Paredes-Fortuny et al. 2015), X-ray (Li et al. 2011, 2012; Chernyakova et al. 2017; D'A'ı et al. 2016), GeV γ -ray (Hadasch et al. 2012; Ackermann et al. 2013; Xing, Wang, & Takata 2017; Chernyakova et al. 2023), and possibly TeV γ -ray (Ahnen et al. 2016). The long-term multiwavelength modulation properties have been investigated in great detail by Chernyakova et al. (2012), Saha et al. (2016), and Jaron (2021). In particular, Chernyakova et al. (2012) found that the X-ray peak shows a systematic drift from the orbital phase of 0.35 to 0.75 within one superorbital period. They also discovered a constant time delay between the X-ray and radio peaks, indicating different emitting regions in these two bands. A detailed study of systematic phase-offsets in radio, X-ray, GeV, and TeV γ -rays are also presented in Jaron (2021). Saha et al. (2016) found that long-term flux variations as seen by Fermi /LAT, RXTE /PCA, Swift /XRT, and Owens Valley Radio Observatory can be well fitted by a sinusoidal function of the superorbital period, with peak positions and amplitudes varying significantly with the binary phase.  \nThe origin of superorbital modulations in LS I+61 · 303 is far from well understood. Zamanov et al. (1999) suggested that long-term modulations could be due to cyclical changes in the Be star's mass-loss rate or disk density. Alternatively, under the microquasar scenario, Massi & Torricelli-Ciamponi (2014) proposed that the superorbital period is caused by the beat frequency between the orbital motion of the CO with P orb = 26 . 496 d and the precession of the microquasar jet with a period of P jet = 26 . 926 d. However, this scenario may be incompatible with the detection of radio pulsations and the lack of accretion signals in X-rays (e.g., Sidoli et al. 2006; Chen & Takata 2022; Saavedra et al. 2023).  \nSuperorbital periods are seen in many X-ray binaries (XRBs; e.g., Sood et al. 2007; Kotze & Charles 2012; Townsend & Charles 2020) but are rare in γ -ray binaries. LS I+61 · 303 is the only γ -ray binary with long-term modulations discovered so far. For Be/XRBs, it has been suggested that the optical superorbital variability is related to the formation and depletion of the Be disk (Alcock et al. 2001; Suffak, Jones, & Carciofi 2022) or the precession and warping of the disk due to the tidal effect of the CO companion (Rajoelimanana, Charles, & Udalski 2011). Those processes can modulate accretion rates onto the CO, leading to long-term modulations of X-rays. Alternatively, Martin (2023) proposed that the stellar spin-axis precession can also drive superorbital modulations in Be/XRBs, particularly for those with shorter orbital periods.  \nLS I+61 · 303 shares similar properties as those Be/XRBs, although the detection of radio pulsation indicates its broadband emission is more likely to be powered by the rotational pulsar instead of the accretion process. In this paper, we will investigate the disk precession scenario for the superorbital modulations of LS I+61 · 303 under the interaction model between the pulsar wind and stellar outflows. We only focus on the X-ray emission since it is widely believed to be produced by synchrotron radiation in the intrabinary shock (IBS), while emission processes in other wavebands (radio and γ -rays) are much more complicated and still underdebated. The paper is organized as follows. Section 2 describes the disk precession geometry and the emission model for LS I+61 · 303. In Section 3, we present expected X-ray flux modulations over the orbital and superorbital periods and compare observational data. We plot the Corbet diagrams of LS I+61 · 303 with other Be/XRBs and discuss possible physical origins that drive the precession of the disk in Section 4. Finally, a brief summary of our emission model is provided in Section 5, and we also discuss the observational properties in other wavebands.", '2.1. Precessing disk': "The stellar outflows of Be stars are characterized by two components: a fast-moving dilute polar wind and a slowmoving dense equatorial disk (Waters 1986). The exact density and velocity distribution profiles of stellar outflows can be very complicated, particularly for the decretion disk (e.g. Torres et al. 2012). Here, we adopt a simplified outflow model for Be stars used in Petropoulou et al. (2018) and Chen et al. (2019). Assuming that stellar outflows are axisymmetric along the disk normal, we model the ram pressure of stellar outflows as:  \np w ( R,θ ) = p 0 R -2 (1 + G | cos θ | m ) , (1)  \nFigure 1. Schematic diagram of the stellar disk precession scenario (not to scale). The normal vector of the stellar disk (solid green arrow) rotates along the precession axis (solid purple arrow) with the trajectory indicated by the dashed green circle.  \n<!-- image -->  \nwhere R is the radial distance from the star, θ is the latitude angle measured from the stellar equator, and p 0 is the polar wind ram pressure. The disk is described by the equator-to-pole pressure contrast G and the confinement index m . A higher value of m means a well-confined disk (Chen & Takata 2022).  \nWe attribute superorbital modulations of LS I+61 · 303 to the interactions of pulsar wind with a precessing Be disk. As illustrated in Figure. 1, we define a Cartesian coordinate system with the origin on the stellar center, x -axis toward the periastron, and z -axis along the orbital axis. Then, the unit vectors of the pulsar and the disk precession axis are (e.g., Chen et al. 2021a):  \ne psr = R [ e z , φ ] e x =(cos φ, sin φ, 0) ,  \ne pre = R [ e z , φ p ] R [ e y , -( π/ 2 -i p )] e x =(sin i p cos φ p , sin i p sin φ p , cos i p ) ,  \n(2) (3)  \nwhere R [ e i , ξ ] is the rotation matrix around the unit vector e i with the counterclockwise rotated angle of ξ , φ is the NS's true anomaly, i p and φ p are the inclination angle and the true anomaly of the precession axis, respectively. The normal vector of the precessing disk at a given time t is  \ne disk ( t ) = R [ e pre , ˙ ω ( t -t 0 )] e disk , 0 , (4)  \nwith  \ne disk , 0 = R [ e z , φ d , 0 ] R [ e y , -( π/ 2 -i d , 0 )] e x =(sin i d , 0 cos φ d , 0 , sin i d , 0 sin φ d , 0 , cos i d , 0 ) , (5)  \nbeing the direction of the disk normal at the reference time t 0 , ˙ ω = 2 π/P pre , and P pre is the precession period.  \n2.2. Intrabinary shock  \nIn γ -ray binaries, the relativistic pulsar wind will be terminated by the intense outflows from the massive star, forming a pair of shocks. The contact discontinuity of shocks is determined by the dynamical balance of two winds:  \nL sd 4 πr 2 s c = p 0 R -2 s (1 + G | cos θ | m ) , (6)  \nwith L sd being the NS's spin-down power, c is the light speed, and the latitude angle is given by:  \nθ ( t ) = π/ 2 -arccos[ e disk ( t ) · e psr ( φ )] . (7)  \nIn the above equation, we need to express the true anomaly φ as a function of time t (Battin 1999):  \nφ ( t ) = M +2 ∞ ∑ k =1 1 k [ ∞ ∑ n = -∞ J n ( -ke ) β | k + n | e ] sin kM, (8)  \nM ( t ) = M 0 + 2 π P orb ( t -t 0 ) , (9)  \nwhere β e = e/ (1 + √ 1 -e 2 ), e is the orbital eccentricity, J n are Bessel functions of the first kind, and M 0 is the mean anomaly at the epoch of t 0 .  \nIt is useful to introduce the momentum ratio of two winds as (Chen et al. 2019):  \nη ( t ) = L sd /c 4 πp 0 (1 + G | cos θ ( t ) | m ) , (10)  \ntherefore, the distances of the IBS apex from the NS and star can be simply written as:  \nr s ( t ) = d ( φ ) η ( t ) 1 / 2 1 + η ( t ) 1 / 2 , (11)  \nand  \nrespectively. The binary separation is  \nwith a being the semi-major axis.", '2.3. Synchrotron radiation': 'The IBS will compress the pulsar wind magnetic field in the downstream region and accelerate the charged particles (which are mainly electron-positron pairs, hereafter e ± pairs) to VHE energies. The accelerated e ± pairs move in the magnetic field and emit synchrotron photons.  \nThe magnetic field strength in the downstream of IBS is:  \nB ( t ) = ( 2 /epsilon1 B L sd r s ( t ) 2 c ) 1 / 2 , (14)  \nwhere the magnetic equipartition parameter /epsilon1 B is related to the pulsar wind magnetization factor σ as /epsilon1 B /similarequal 9 σ/ 2 for σ < 0 . 1 or /epsilon1 B /similarequal 1 / 2 for σ > 0 . 1 (Kennel & Coroniti 1984). It is obvious that a smaller value of r s during the disk passages will lead to an increase in the magnetic field of IBS.  \nThe e ± pairs from pulsar wind are expected to be accelerated into a power-law distribution. We write the pair injection rate as (Kirk, Ball, & Skjæraasen 1999):  \nQ ( γ ) /similarequal /epsilon1 e L sd ( p -2) γ 2 -p min m e c 2 · γ -p , γ ≥ γ min , (15)  \nR s ( t ) = d ( φ ) 1 1 + η ( t ) 1 / 2 , (12)  \nd ( φ ) = a 1 -e 2 1 + e cos φ , (13)  \nwhere /epsilon1 e is the fraction of the pulsar luminosity being transferred to accelerated e ± pairs, and m e is electron mass. Assuming that the unshocked pulsar wind contains mono-energetic e ± pairs with Lorentz factor of γ w , the minimum Lorentz factor of accelerated e ± pairs is γ min = γ w ( p -2) / ( p -1). Since we only consider the expected emission in soft X-ray band (i.e., 0 . 3 -10 keV), emitting particles considered here are dominated by adiabatic cooling with the timescale of (Dubus 2006; Kong, Cheng, & Huang 2012)  \nτ ad ( t ) /similarequal 3 r s ( t ) c , (16)  \nand the evolved distribution of cooled e ± pairs can be simply written as (Kirk, Ball, & Skjæraasen 1999):  \nN ( γ ) /similarequal Q ( γ ) τ ad . (17)  \nUnder the monochromatic approximation, the synchrotron emissivity of a single electron is (Petropoulou et al. 2018):  \nP ν ( γ ) /similarequal 4 3 cσ T u B γ 2 δ ( ν -ν c ) , (18)  \nwhere  \nhν c /similarequal 3 2 B B cr γ 2 m e c 2 , (19)  \nis the characteristic synchrotron energy, u B = B 2 / 8 π , B cr = 4 . 414 × 10 13 Gs is the quantum critical field strength, and σ T is the Thompson cross section. The specific synchrotron luminosity is obtained by integrating over the particle distribution:  \nL  \nν  \n=  \n∫  \nN  \n(  \nγ  \n)  \nP  \nν  \n(  \nγ  \n)  \ndγ  \nwhich gives the analytical expression of synchrotron luminosity in the X-ray frequency range of [ ν X1 , ν X2 ] as:  \nL X ( t ) /similarequal C f X /epsilon1 e /epsilon1 1+ p 4 B L 5+ p 4 sd γ p -2 w r s ( t ) 1 -p 2 , (20)  \nwhere C = 2 17 -3 p 4 3 p -5 2 π 3 -p 2 c -19 -3 p 4 m -p -3 2 e q p +5 2 e ( p -2) p -1 ( p -1) 2 -p , f X = ( ν 1 -β X1 -ν 1 -β X2 ) / ( β -1), and β = ( p -1) / 2 (Petropoulou et al. 2018). Obviously, the modulation of X-ray luminosity is completely determined by the variation of the shock distance from the NS. Due to orbital motion and disk precession, interactions between pulsar wind and stellar outflows lead to cyclic changes in r s ( t ), resulting in the orbital and superorbital modulations of X-ray flux.', '3. RESULT': "In this section, we will use the emission model described above to calculate the time-dependent of X-ray flux, and refold it into the orbital and superorbital phases as:  \nΦ orb = t -t 0 P orb -int ( t -t 0 P orb ) , (21)  \nΘ sup = t -t 0 P sup -int ( t -t 0 P sup ) , (22)  \nwhere int( x ) is the intergerpart of x . According to Eq.(20), the synchrotron luminosity depends on the following parameters: /epsilon1 e , /epsilon1 B , L sd , γ w , p and r s . Among them, only r s is time-dependent, which is determined by the orbital solution and the outflow properties of binary stars. Unfortunately, the orbital elements of LS I+61 · 303 are poorly determined, and several orbital solutions are obtained (see Fig. 2). The observed linear drift of X-ray peak from orbital phase 0.35 to 0.75 in one superorbital period (Chernyakova et al. 2012) favors a less elliptical orbit with the periastron phase falling within this range. Therefore, in calculations, we will use the orbital parameters given by Kravtsov et al. (2020), which specify an eccentricity of e /similarequal 0 . 06 and a periastron phase of Φ p ∼ 0 . 6. These values were derived from fitting the star's optical linear polarization, although they differ significantly from those obtained  \nFigure 2. Orbits of LS I+61 · 303 with orbital solutions of Hutchings & Crampton (1981); Casares et al. (2005); Grundstrom et al. (2007); Aragona et al. (2009); Kravtsov et al. (2020) as labeled in the plot. The sizes of the star and orbits are to scale.  \n<!-- image -->  \nthrough measurements of radial velocity variations in spectral lines. For stellar outflows, we take the typical value of polar wind ram pressure as p 0 = 2 × 10 26 g · cm · s -2 , the pressure contrast and confinement parameter of the disk as G = 100 and m = 100, respectively. The effects of these disk parameters on light curves have been explored in Chen et al. (2019). The other model parameters we adopt are as follows: the NS's luminosity L sd /similarequal 5 × 10 36 erg · s -1 , the Lorentz factor of e ± pairs in pulsar wind γ w = 10 4 , and the power-law index of accelerated particles p = 2 . 1. The magnetic equipartition parameter and the fraction of spin-down luminosity transferred by shock-accelerated e ± pairs are /epsilon1 B = 0 . 005 and /epsilon1 e = 0 . 1, respectively. It should be noted that these parameters can be scaled accordingly for other values according to Eq.(20). The disk precession timescale adopted in calculations is P pre = 2 P sup because the emission is expected to be the same when the disk is at the longitude of ascending node Ω and at Ω + 180 · (e.g. Martin 2023). The observational and model parameters are listed in Table. 1.  \nAs the NS orbits around the Be companion, the relativistic pulsar wind interacts with the inclined disk twice per orbit. The additional disk pressure pushes the IBS closer to the NS, enhancing the magnetic field and, therefore, producing the asymmetric two-peak profile in the orbital light curve. Meanwhile, the disk precession leads to cyclical changes in the magnetic field of IBS at each orbital phase, which results in a long-term flux modulation at given orbital phases. The most uncertain parameters are the inclination angles of the disk normal and the precession axis, which may strongly affect the shapes of light curves. In Fig. 3, we explore how these angles affect the orbital and superorbital modulations of X-ray flux at the phases of Θ sup = 0 . 5 (top panel) and Φ orb = 0 . 5 (bottom panel), respectively. Obviously, in the case of i p = 0 · and i d , 0 = 0 · , the disk is coplanar with the orbital plane with no precession, then the orbital light curve simply peaks in periastron without long-term modulation at any given orbital  \nFigure 3. Calculated X-ray flux over the orbital period (top panel) and superorbital period (bottom panel) with different values of i p and i d , 0 at the phases of Θ sup = 0 . 5 and Φ orb = 0 . 5, respectively. For clarity, we do not include the Doppler boosting effect in this plot.  \n<!-- image -->  \nFigure 4. Contour plot of synchrotron X-ray flux as a function of orbital and superorbital phases. The horizontal and vertical histograms show the calculated superorbital light curve at the phase of Φ orb = 1 . 0 and the calculated orbital light curve at the phase of Φ sup = 1 . 0, respectively. The black dashed lines are the observed orbital phase drifting trend of X-ray peak along the superorbital period given by Chernyakova et al. (2012).  \n<!-- image -->  \nphase. With increases of i p or i d , 0 , the asymmetric double-peak features in orbital light curves arise (top panel). For superorbital light curves (bottom panel), higher inclination angles of the disk normal or the precession axis mean that the NS interacts with the precessing disk only at certain epochs, and the X-ray flux at a given orbital phase remains constant during most of the rest superorbital period. Furthermore, a higher misaligned angle between the precession and orbital axes would also lead to a plateau structure or two sinusoidal peaks in the superorbital light curves.  \nThe X-ray superorbital light curves of LS I+61 · 303 show a very smooth sinusoidal-like modulation (e.g., Li et al. 2012), suggesting that the Be disk is tilted to the orbital plane with a relatively small angle ( i d ∼ 10 · ), and the disk normal vector undergoes precession around the orbital axis ( i p ∼ 0 · ). This can be explained by the fact that if a disk is misaligned, the CO's tidal torque will tend to align it to the orbital plane (e.g., Okazaki 2017). In Fig.4, we show the contour plot of synchrotron-powered X-ray flux as a function of orbital and superorbital phases, with the disk parameters listed in Table. 1. The horizontal and vertical plots are the corresponding superorbital and orbital light curves at the phases of Φ orb = 1 . 0 and Θ sup = 1 . 0, respectively, as indicated in the main diagram. The black dashed lines are the observed orbital phase drifting trend of the X-ray peak along the superorbital period given by Chernyakova et al. (2012). The corresponding calculated X-ray light curves over orbital and superorbital periods at each phase with comparisons of Swift /XRT data are presented in Fig. 5. In calculations, we also include the Doppler boosting effect of shocked flows (e.g. Dubus, Cerutti, & Henri 2010) with a very mildly-moving velocity ( v sh /similarequal 0 . 04 c ) to give better fittings of observations. The overall phase drifting trend and calculated light curves are roughly consistent with the observations, although there are still some discrepancies (i.e., the superorbital light curves at the phases of  \nFigure 5. Calculated X-ray flux over orbital period (top panel) and superorbital period (bottom panel) at different phases with comparisons of Swift /XRT data taken from Saha et al. (2016).  \n<!-- image -->  \nΦ orb = 0 . 5 -0 . 7). The discrepancies could be due to our oversimplifications of the stellar outflow. In particular, the rotation of the disk can induce additional changes of shock geometry and, therefore, to the Doppler boosting flux when the NS is moving around inferior conjunction (i.e., Φ orb ∼ 0 . 6 -0 . 7). Furthermore, photoelectric absorption of soft X-rays and inhomogeneities in the stellar outflows can also affect the observed flux (Szostek & Dubus 2011; Chernyakova et al. 2017).", '4. SUPERORBITAL PERIODS IN BE/X-RAY BINARIES': "Long-term flux modulations are seen in many types of XRBs. In Fig. 6, we plot the Corbet diagrams of LS I+61 · 303 and other Be/XRBs with measured values of P NS , P orb and P sup , which are listed in Table. 2. In the P NS -P orb diagram (top panel), the NS spin-periods of Be/XRBs follow a loose relationship of P NS ∝ P 2 orb as already noted by Corbet (1986). This relationship can be explained by the NS's equilibrium spin-period with a torque of accreted material from the Be disk exerted on the NS (Waters & van Kerkwijk 1989):  \nP eq /similarequal 2 3 / 2 πµ 6 / 7 δ -9 / 7 (2 GM X ) 11 / 7 ( πρ 0 v 3 0 ) -3 / 7 ( GMP 2 orb 4 π 2 R 2 /star ) (4 n -6) / 7 (23)  \nwhere µ is the NS's magnetic moment; δ = v d , r /v rel , v d , r is the disk radial velocity, v rel is the disk velocity relative to the NS's orbital motion; M = M X + M /star is the total binary mass, ρ 0 and v 0 are the base density and velocity of the disk at R /star , respectively (Li & van den Heuvel 1996). The observed slope in the P NS -P orb diagram around 2 indicates the radial density distribution index of the disk being around n /similarequal 3 . 25. For LS I+61 · 303, detections of radio pulsations (Weng et al. 2022) and lack of accretion signatures in X-rays (Sidoli et al. 2006) suggest that the NS is in the rotation-powered state, and the equilibrium relationship may not hold for the system. This may explain why LS I+61 · 303 stands as an outlier in the P NS -P orb diagram.  \nThe superorbital periods in Be/XRBs have been suggested to be related to the formation and depletion of Be disks (Alcock et al. 2001). However, the linear correlation in the P sup -P orb diagram suggests a binary dynamical origin of the superorbital modulations. For Be/XRBs, the evolution of the Be disk would be significantly affected by its companion object. The torque on the disk from the NS companion is (e.g. Lai 2014):  \nT d /similarequal GM X M d R 2 d 4 a 3 cos i d , (24)  \nand the precession rate of the disk normal around the binary axis can be estimated as:  \n˙ ω d /similarequal 1 4 M X M /star ( GM /star a 3 ) 1 / 2 ( R d a ) 3 / 2 cos i d , (25)  \nwhere M d and R d are the mass and size of the disk. Assuming that R d has the same order of the binary separation, and taking the parameters of LS I+61 · 303 into the above equation, we can estimate the disk precession timescale with P d ∼ 1000 d, which agrees with LS I+61 · 303 within a factor of 3. In addition to the disk precession, the NS companion also exerts a torque on the oblate Be star, resulting in the stellar spin-axis precession (Kuhnel et al. 2017; Martin 2023). In this case, the disk can also precess like a solid body if its viscosity parameter is high enough. However, this precession timescale is more than 100 yr for LS I+61 · 303, which can be neglected.", '5. DISCUSSION AND CONCLUSION': 'LS I+61 · 303 is a unique source for us to investigate dynamical interactions within binary systems and the outflow properties of massive stars. The most distinctive feature of this binary is its superorbital modulations seen in nearly all wavelengths. The detection of radio pulsations and the lack of accretion signals strongly favor the wind-interaction scenario instead of the microquasar model. Assuming synchrotron radiations of the IBS produce X-rays, we calculate the expected flux modulations folded on the orbital and superorbital periods with comparisons of Swift/XRT data. The asymmetric double-peak profiles, as seen in the orbital light curves of X-rays, are attributed to the interaction between the pulsar wind and the slightly inclined Be disk. Considering the disk precession, the sinusoidal-like long-term modulations in X-rays are reproduced. The linear correlation of LS I+61 · 303 with other Be/XRBs in the P sup -P orb diagram suggests that the Be disk precession is induced by the tidal effects of the NS companion.  \nIn this paper, we have only focused on the soft X-ray emission properties detected by Swift/XRT, whereas superorbital modulation has been seen in nearly all wavelengths. In particular, Chernyakova et al. (2012) performed a  \nFigure 6. Corbet diagrams (top panel: P NS -P orb ; bottom panel: P sup -P orb ) of LS I+61 · 303 and other Be/XRBs with superorbital periods as listed in Table. (2). The black lines in each panel represent the empirical fitting functions with P NS ∝ P 2 orb and P sup ∝ P orb , respectively.  \n<!-- image -->  \nlong-term multiwavelength analysis (including radio, hard X-ray, and GeV γ -ray) of LS I+61 · 303, and discovered a constant time lag between the radio and X-ray peak. It implies different emitting mechanisms and/or regions in these two bands. It has been suggested that radio emission is likely to be produced by cooled particles at the shock tail at a much larger distance (Mold\'on, Rib\'o, & Paredes 2013), and the shock could be deflected due to the orbital motion and disk rotation. Therefore, the radio emission is directly related to the large-scale shock structure, fluid propagation, and cooling processes of particles. Additionally, absorption processes within the binary (including free-free absorption and synchrotron self-absorption) will also affect the observed radio flux (Chen et al. 2021a). These factors may explain why the system shows a different modulation pattern compared to X-rays.  \nLS I+61 · 303 also shows energy-dependent variability in the GeV γ -rays according to the analysis of over 14 years of Fermi/LAT data (Chernyakova et al. 2023). In particular, the low-energy (0.1-0.3 GeV) flux is strongly modulated with the orbital period, while modulation in the high-energy (1-10 GeV) is strongly suppressed. It has been suggested that GeV emissions from gamma-ray binaries are likely contributed by several components, including outer-gap emission in pulsar magnetosphere, anisotropic IC scattering within cold pulsar wind and shock emission (e.g., Takata et al. 2014; Chen et al. 2021b). The emissions from the latter two can be significantly affected by the disk precession. To understand the variability pattern of GeV emission, all of those emission components should be taken into account properly.  \nSince its discovery over 45 years ago, LS I+61 · 303 has been the subject of intensive observational campaigns across multiple telescopes, resulting in a substantial accumulation of data. This mysterious system exhibits numerous unexpected characteristics that set it apart from other γ -ray binaries. Despite extensive research efforts, a comprehensive theoretical framework capable of elucidating the multiwavelength flux modulation - spanning from radio to TeV energies - remains elusive. Although our current model provides a tentative approach to explain the long-term flux variations in binaries, a more sophisticated emission model should be developed to explain the multiwavelength data.  \n- We are very grateful to the referee for insightful comments on the manuscript and to Dr. Lab Saha for sharing the data 1\n- of LS I+61 · 303. We also thank Prof. Dong Lai and Prof. Weihua Lei for helpful discussions. Y.W.Y. is supported 2\n- by the National Key R&D Program of China (No. 2021YFA0718500), the National SKA Program of China (No. 3\n- 2020SKA0120300), the China Manned Spaced Project (No. CMS-CSST-2021-A12) and the National Natural Science 4\n- Foundation of China (No. 12393811). J.T. is supported by the National Key Research and Development Program of 5\n- China (No. 2020YFC2201400) and the National Natural Science Foundation of China (Nos. U1838102, 12173014). 6\n- A.M.C. is supported by the China Postdoctoral Science Foundation (No. 2023T160410). 7  \nTable 1. Parameters of LS I+61 · 303.  \nReference. (1) Kravtsov et al. 2020; (2) Gregory 2002; (3) Massi & Torricelli-Ciamponi 2016.  \nTable 2. List of Be/XRBs with superorbital modulations.  \nReference. (1) Weng et al. 2022; (2) Gregory 2002; (3) Massi & Torricelli-Ciamponi 2016; (4) Doroshenko et al. 2021; (5) Richardson et al. 2023; (6) Corbet & Krimm 2013; (7) Islam et al. 2015; (8) Doroshenko, Tsygankov, & Santangelo 2018; (9) Liu et al. 2022 ; (10) White et al. 1980; (11) Alfonso-Garz\'on et al. 2017; (12) Vasilopoulos et al. 2016; (13) Vasilopoulos et al. 2014; (14) Rajoelimanana et al. 2017; (15) Rajoelimanana, Charles, & Udalski 2011.  \n- D\'A\'ı, A., Cusumano, G., La Parola, V., Segreto, A., &  \nMineo, T. 2016, MNRAS, 456, 1955. Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, VI Microquasar Workshop: Microquasars and Beyond, 52.1. Doroshenko, V., Tsygankov, S., & Santangelo, A. 2018, A&A, 613, A19. Doroshenko, V., Santangelo, A., Tsygankov, S. S., & Ji, L. 2021, A&A, 647, A165. Dubus, G. 2006, A&A, 456, 801. Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18. Gregory, P. C., & Taylor, A. R. 1978, Nature, 272, 704. Gregory, P. C. 2002, ApJ, 575, 427. Gregory, P. C., & Neish, C. 2002, ApJ, 580, 1133. Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437. Gupta, S., & Bottcher, M. 2006, ApJL, 650, L123. Hadasch, D., Torres, D. F., Tanaka, T., et al. 2012, ApJ, 749, 54. Hutchings, J. B., & Crampton, D. 1981, PASP, 93, 486. Islam, N., Maitra, C., Pradhan, P., & Paul, B. 2015, MNRAS, 446, 4148. Jaron, F., Kiehlmann, S., & Readhead, A. C. S. 2024, A&A, 683, A228. Jaron, F., & Massi, M. 2014, A&A, 572, A105. Jaron, F. 2021, Universe, 7, 245. Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694. Kirk, J. G., Ball, L., & Skjæraasen, O. 1999, Astroparticle Physics, 10, 31. Kravtsov, V., Berdyugin, A. V., Piirola, V., et al. 2020, A&A, 643, A170. Kong, S. W., Cheng, K. S., & Huang, Y. F. 2012, ApJ, 753, 127. Kotze, M. M., & Charles, P. A. 2012, MNRAS, 420, 1575. Kuhnel, M., Rothschild, R. E., Okazaki, A. T., et al. 2017, MNRAS, 471, 1553. Lai, D. 2014, MNRAS, 440, 3532. Li, J., Torres, D. F., Zhang, S., et al. 2011, ApJ, 733, 89. Li, J., Torres, D. F., Zhang, S., et al. 2012, ApJL, 744, L13. Li, X.-D., & van den Heuvel, E. P. J. 1996, A&A, 314, L13. Lipunov, V. M., & Nazin, S. N. 1994, A&A, 289, 822. Liu, W., Yan, J., Reig, P., et al. 2022, A&A, 666, A110. Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1P. Marcote, B., Rib\'o, M., Paredes, J. M., et al. 2016, MNRAS, 456, 1791. Marcote, B. 2017, XII Multifrequency Behaviour of High Energy Cosmic Sources Workshop (MULTIF2017), 45. Martin, R. G. 2023, MNRAS, 523, L75. Massi, M., & Torricelli-Ciamponi, G. 2014, A&A, 564, A23.  \nThis figure "orcid-ID.png" is available in "png"GLYPH<10> format from:  \nhttp://arxiv.org/ps/2409.04818v1'}