The Puzzling Superorbital Period Variation of the Low-mass X-ray Binary 4U 1820-30

In the power spectral analysis, all the light curves were rebinned into daily averages for consistency. All the data points with a signal-to-noise ratio less than 2$\sigma$ were filtered out for further analysis. To probe the superorbital periods of various observations, the Lomb-Scargle (LS) periodogram (Scargle, 1982) was applied to generate the power spectra. The errors of signal frequencies were estimated by the method proposed by Horne & Baliunas (1986):

The power spectra are depicted in Figure 2, and the detected superorbital periods are outlined in Table 1. It is apparent that the superorbital period deviates from the expected stability suggested by triple model, showingt a tendency to decrease over time. By incorporating the superorbital period reported by (Priedhorsky & Terrell, 1984) from the Vela5B light curve, we estimated the timescale of the superorbital change by fitting a linear function to the detected superorbital periods over time (see Figure 3). This result in a period derivative of $\dot{P}_{sup}/P_{sup}=(-7.37\pm 0.33)\times 10^{-4}$ $yr^{-1}$, corresponding to an evolution timescale of 1,357 years. This period derivative exceeds the upper limit proposed by CG10 ($|\dot{P}_{sup}/P_{sup}|<2.2\times 10^{-4}$ $yr^{-1}$). However, the linear fitting yielded a reduced $\chi^{2}$ of 10.2, suggesting that the superorbital evolution of 4U 1820-30 is likely more complex than the constant period derivative model suggests. Further variations in the superorbital period from phase analysis will be demonstrated in Section 3.2.

Additionally, to compare the amplitudes of superorbital modulation, we folded these five light curves using two kinds of linear ephemerides. The first one is the optimal ephemeris proposed by CG01 (hereafter CG01 ephemeris),

In this research, we aimed to trace the long-term evolution of the superorbital phase of 4U 1820-30, necessitating the analysis of superorbital phases measured from different instruments. However, time lags between different energy bands are often observed in astronomical time series. For instance, soft phase lags of pulsed emissions are commonly noted in accreting millisecond X-ray pulsars (Cui et al., 1998; Patruno & Watts, 2021). Hence,a coherence test was conducted to verify if there was a significant time lag between any of two light curves from different instruments. However, this test could be only performed on the light curves with overlapping observation times. Among the five light curves we analyzed in this study, except for Ginga ASM, other 4 light curves had overlapping observation times for each other, resulting in 6 pairs of light curves for the coherence test. For each pair of light curves, only overlapping parts were selected for coherence test. The power spectra were obtained the superorbital periods for the corresponding light curves. The superorbital modulation profiles of both light curves were conducted by folding the mean period measured from the power spectra with an arbitrary but fixed phase zero epoch. We discovered that all the profiles could be well fitted with a four-component sinusoidal function, that is, $r(\phi)=a_{0}+\sum_{k=1}^{4}[a_{k}cos(2\pi k\phi)+b_{k}sin(2\pi k\phi)]$. To measure the possible time delay between the two instruments, we applied the cross-correlation between the best-fitted modulation profiles. The test results are shown in Table 2. The phase difference is generally no more than 0.026 cycle, which is much smaller the phase jitters ($\sim$0.1 cycle, see CG01). We conclude that there is no significant time delay is observable among these 4 instruments for the superorbital modulation of 4U 1820-30.

To trace the evolution of superorbital phase, we segmented the light curves and folded them to derive the modulation profiles. For instruments highly sensitive to superorbital modulation, like RXTE ASM and MAXI light curves, two cycles ($2\times 171$ days) per segment sufficed to yield clear profiles. In the case of Swift BAT observations, where no superorbital modulation was detected for the photon energy larger than 24 keV (Farrell et al., 2009), we selected four cycles as a data segment to ensure significant profile detection. As for the RXTE PCA light curve, due to the observation gaps, we adopted four cycles per segment to create the profiles. However, only three data segments provided sufficient phase coverage for further analysis. Given the very low sensitivity of Ginga ASM, a clear profile could only be obtained by folding the entire light curve.

Following the approach of CG01, we selected the superorbital minimum as the fiducial point of the superorbital phase. Ideally , we would fold a light curve segment using a fix ephemeris, such as CG01 ephemeris (Eq. 3) to determine the phase (i.e. $\phi_{CG01}$). However, as indicated in Section 3.1, the CG01 ephemeris is unlikely to be an optimal ephemeris for the all observations, particularly for recent ones (e.g. Swift BAT and MAXI observations), which could lead to profile deformation. To precisely determine the $\phi_{CG01}$, we folded the light curve segments using the optimal linear ephemeris specific to each instrument (local ephemeris). This involved folding the data by the period obtained from power spectral analysis (see Table 1) and an arbitrary but fixed phase zero epoch. A typical modulation profile is depicted in Figure 4. The phase ($\phi_{local}$) of a data segment was determined by fitting a four-component sinusoidal function and identifying the phase corresponding to the intensity minimum (fiducial point). This phase value, along with the local ephemeris, facilitated the evaluation the superorbital minimum time ($t_{m}$) closest to the mid of observation time of the data segment. Subsequently, $t_{m}$ was then folded by CG01 ephemeris (Eq. 3) to obtain the phase $\phi_{CG01}$.

The superorbital orbital phases ($\phi_{CG01}$) are listed in Table 3, and their evolution is illustrated in Figure 5. It is noteworthy that while the superorbital modulation displays strong periodicity in the power spectra, the modulation profile varies from cycle to cycle, exhibiting the quasi-periodic nature as described in  Zdziarski et al. (2007a). This variability induces phase jitters of $\sim$0.1 cycle, evident in Figure 5 and CG01. These phase jitters are considerably larger than the error estimated from photon statistics ($\sim$0.005 cycle). Despite the presence of phase jitters, a discernible phase evolution trend can be discerned in Figure 5. However, independently evaluating errors from phase jitters is difficult, which depend on the evolution model. In the subsequent analysis, we utilized the unweighted fitting method outlined by Press et al. (2002) to update the ephemeris for the superorbital modulation of 4U 1820-30.

As previously mentioned in Section 1, Type-I X-ray bursts of 4U 1820-30 are exclusively observable during the low state (Clark et al., 1977; Stella et al., 1984). Further confirmed by CG01, indicated the Type-I X-ray bursts were detected only within $\pm$23 days around superorbital minima for bursts reported before 1985. This supports the notation that the superorbital modulation stems from changes in the accretion rate changes rather than occultation effects. However, CG01’s statistics only included four burst active dates. Subsequently, more X-ray bursts of 4U 1820-30 were detected. With the updated superorbital ephemerides, this evidence can be further substantiated. Although the possibility exists that the X-ray bursts occur in another low state, which may deviate significantly from the superorbital minima (e.g.the brief low state found by  Šimon (2003)), it is likely that most of X-ray burst active times would cluster around the superorbital minima.

In this study, we collected the reported burst active dates of 4U 1820-30 from Grindlay et al. (1976); Vacca et al. (1986); Haberl et al. (1987); Zdziarski et al. (2007b); García et al. (2013), and (Yu et al., 2024), as well as four superbursts discovered by Strohmayer & Brown (2002); in’t Zand et al. (2011); Serino et al. (2021a), and (Serino et al., 2021b). The deviations in burst active dates relative to the nearest superorbital minima predicted by CG01 and three ephemerides derived in this study are shown in Figure 8. It is evident that large deviations can be observed after the low power state (MJD 52630) when using the superorbital minima predicted by the CG01 ephemeris. The root-mean-square deviations are 41.2, 32.7 27.4 and 27.0 days for CG01, linear, glitch and quadratic ephemerides, respectively. This provides supportive evidence that glitch and quadratic ephemerides are better than CG01 and linear ephemerides in describing the superorbital phase evolution of 4U 1820-30. Notably, X-ray bursts occurring between MJD 50399 and MJD 50343 exhibited large deviations ($\sim 60$ days) for all ephemerides. Upon examining the RXTE ASM light curve, we observed that these X-ray bursts occurred around a BLS with a count rate of only 12.0 cts/s compared to the mean count rate of 21.0 cts/s. This observation reaffirms that superorbital modulation is primarily caused by accretion variations rather than occultation effects.

As previously mentioned in Section 1, Type-I X-ray bursts of 4U 1820-30 are exclusively observable during the low state (Clark et al., 1977; Stella et al., 1984). Further confirmed by CG01, indicated the Type-I X-ray bursts were detected only within $\pm$23 days around superorbital minima for bursts reported before 1985. This supports the notation that the superorbital modulation stems from changes in the accretion rate changes rather than occultation effects. However, CG01’s statistics only included four burst active dates. Subsequently, more X-ray bursts of 4U 1820-30 were detected. With the updated superorbital ephemerides, this evidence can be further substantiated. Although the possibility exists that the X-ray bursts occur in another low state, which may deviate significantly from the superorbital minima (e.g.the brief low state found by  Šimon (2003)), it is likely that most of X-ray burst active times would cluster around the superorbital minima.

In this study, we collected the reported burst active dates of 4U 1820-30 from Grindlay et al. (1976); Vacca et al. (1986); Haberl et al. (1987); Zdziarski et al. (2007b); García et al. (2013), and (Yu et al., 2024), as well as four superbursts discovered by Strohmayer & Brown (2002); in’t Zand et al. (2011); Serino et al. (2021a), and (Serino et al., 2021b). The deviations in burst active dates relative to the nearest superorbital minima predicted by CG01 and three ephemerides derived in this study are shown in Figure 8. It is evident that large deviations can be observed after the low power state (MJD 52630) when using the superorbital minima predicted by the CG01 ephemeris. The root-mean-square deviations are 41.2, 32.7 27.4 and 27.0 days for CG01, linear, glitch and quadratic ephemerides, respectively. This provides supportive evidence that glitch and quadratic ephemerides are better than CG01 and linear ephemerides in describing the superorbital phase evolution of 4U 1820-30. Notably, X-ray bursts occurring between MJD 50399 and MJD 50343 exhibited large deviations ($\sim 60$ days) for all ephemerides. Upon examining the RXTE ASM light curve, we observed that these X-ray bursts occurred around a BLS with a count rate of only 12.0 cts/s compared to the mean count rate of 21.0 cts/s. This observation reaffirms that superorbital modulation is primarily caused by accretion variations rather than occultation effects.

On the contrary, out of the four superbursts detected, only the one on MJD 51430 (Strohmayer & Brown, 2002) fell within $\pm$23 days region, whereas the other three occurred outside of this timeframe for both glitch and quadratic ephemerdes (see Figure 8). Upon examining the RXTE ASM and MAXI light curves, we observed that the count rates on the dates of superbursts were 57%, 150%, 121% and 138% of the corresponding mean count rates of the light curves for the superbursts detected on MJD 51430 (Strohmayer & Brown, 2002), MJD 55272 (in’t Zand et al., 2011), MJD 59449 (Serino et al., 2021a) and MJD 59543 (Serino et al., 2021b), respectively. It is likely that the low state constraint for the regular X-ray bursts of 4U 1820-30 (Clark et al., 1977; Stella et al., 1984) does not apply to the superbursts. More observations are required for further confirmation.

The initial aim of this study was to further validate the stability of superorbital period of 4U 1820-30, a crucial piece of evidence for the triple model as described in Section 1, which explains its superorbital modulation with a period of $\sim$170 days. Given that 4U 1820-30 resides in NGC 6624, a star-crowded region, previous studies suggested a high likelihood of the binary capturing a third star and forming a stable hierarchical triple system (Grindlay, 1988). Black (1982) proposed a stability criteria for such a triple system as

However, upon analyzing X-ray light curves collected by the sky monitoring/scanning instruments from 1987 to 2023, we discovered a significant change in the superorbital period from 171 days to 167 days, identified through both power spectral analysis (Section 3.1) and phase analysis (Section 3.2) over a time span of $\sim$36 years. This suggests that the ephemeris proposed by CG01 is no longer suitable for describing the superorbital modulation of 4U 1820-30, and the period is not as stable as anticipated by triple model. By analyzing the superorbital phase evolution, we suggested that the superorbital period may have experienced an abrupt change during late 2000 to early 2003 ($T_{g}=$ MJD $52264\pm 466$) or may be constantly changing with a period derivative of $\dot{P}=(-3.58\pm 0.72)\times 10^{-4}$ day/day.

The significant difference between the period detected from Vela5B observation, $176.4\pm 1.3$ days (Priedhorsky & Terrell, 1984), and Ginga observation, $171.12\pm 1.99$ days (see Table 1) suggests that the superorbital period may have experienced another glitch between 1976 and 1987. If 4U 1820-30 is a hierarchical triple system, from Eq. 1, a glitch in superorbital period may be induced by changes in either the binary orbital period or the third star orbital period as

The low power state is a particular phase during the superorbital modulation evolution of 4U 1820-30. The dynamic power spectrum demonstrated in Figure 24 of Kotze & Charles (2012) indicates that in addition to the weaker power detected in the superorbital period, the powers of its first harmonics and a signal of period $\sim$65 days became significant. This suggests that the modulation were more complicate than usual. Our phase analysis results also show that the superorbital phases had larger fluctuation during the low power state (see Figure 5) with an RMSD of 0.11, evaluated by the best glitch model fitting, compared to 0.075 for the phases outside low power state. It is probable that the superorbital period was 171 days before the low power state, became unstable during the low power state, and stabilized at 167 days later. However, it is unclear what the cause of this phenomenon is.

Conversely, the superorbital phase evolution is also well-fitted by a quadratic model, although there is a significant difference (2.2$\sigma$) between the period evaluated as ephemeris extrapolated to the midpoint of Vela5B observation time ($172.84\pm 0.94$ days) and the observed period ($176.4\pm 1.3$ days) as reported by Priedhorsky & Terrell (1984). If 4U 1820-30 is a hierarchical triple system, from Eq. 1, the relation of period derivatives can be written as

If 4U 1820-30 is not a triple system, the constraints regarding the observed value of the binary orbital period derivative may be relaxed. To explain the discrepancy between the positive theoretical value ($\dot{P}_{orb}/P_{orb}>8.8\times 10^{-8}$ yr^-1, Rappaport et al., 1987) and the negative observed value ($\dot{P}_{orb}/P_{orb}=(-5.21\pm 0.13)\times 10^{-8}$ yr^-1, Chou & Jhang, 2023) of the binary orbital period derivative, it has been proposed that 4U 1820-30 is being accelerated by the gravitational potential within the globular cluster NGC 6624 (Tan et al., 1991; Peuten et al., 2014). However, Peuten et al. (2014) suggested that the maximum radial acceleration from the gravitational potential from NGC 6624 itself ($|a_{c,max}/c|=1.3\times 10^{-9}$ yr^-1 is an order of magnitude smaller than the value required to explain the observed binary period period derivative. Therefore, Peuten et al. (2014) proposed three possible scenarios to provide additional acceleration for 4U 1820-30, a flyby stellar mass dark remnant, a intermediate-mass black hole at the center of NGC 6624, and a central concentration of dark remnants. Only the last scenario was preferred because the first two scenarios tend to destroy the triple system (Peuten et al., 2014). However, if 4U 1820-30 is a pure binary system, the first two scenarios become viable explanations for the observed binary orbital period derivative.

The fact that the type-I X-ray bursts of 4U 1820-30 can only be observed in the low state (Clark et al., 1977; Stella et al., 1984) has been reconfirmed in this work (see Section 3.4). This implies that the superorbital modulation of 4U 1820-30 is caused by variations in accretion rate rather than by external absorption or precession of accretion disk. Kotze & Charles (2012) listed eight possible mechanisms to account for the superorbital modulations observed in X-ray binaries, but only the third body (i.e., triple model) and the X-ray state changes can possibly be responsible for the superorbital modulation of 4U 1820-30. If the triple model is ruled out due to the instability of superorbital period, the only remaining mechanism is the X-ray state changes. X-ray state changes refer to variations in mass accretion rate between high and low states due to thermal disk instability, as observed in dwarf novae and soft X-ray transients. Priedhorsky & Terrell (1984) proposed that this mechanism could explain for the superorbital modulation of 4U 1820-30. However, Menou et al. (2002) pointed out that if thermal disk instability could occur in 4U 1820-30 system, the mass transfer rate ${\dot{m}}\leq{\dot{m}_{crit}=4.4\times 10^{16}}$ g s^-1. From the mean flux of $<F_{bol}>=8.7\times 10^{-9}$ erg cm^-2 s^-1 for 4U 1820-30 (Zdziarski et al., 2007b) and the distance of 8.019 kpc for NGC 6624 (Baumgardt & Vasiliev, 2021), we obtained a mean luminosity of $<L>=6.7\times 10^{37}$ erg s^-1 and a mass accretion rate of $\dot{m}_{1}=3.6\times 10^{17}$ g s^-1 for a neutron star with a mass of 1.4 $M_{\sun}$ and a radius of $10^{6}$ cm. It is approximately an order of magnitude larger than the $\dot{m}_{crit}$. Therefore, this mechanism is unlikely to explain the superorbital modulation of 4U 1820-30.

Zdziarski et al. (2007b) discovered that the binary orbital modulation amplitude and the offset phase in 4U 1820-30 depend significantly on the accretion rate, which is highly related to the superorbital modulation phase. The orbital modulation in the X-ray band is believed to be caused by absorption from structures in the disk rim where the accretion flow from the companion impacts the outer edge of the disk (Stella et al., 1987). As the mass loss rate changes, variations of the accretion stream induce changes in the absorption of outer edge structures and the position of impact point. This makes the amplitude and phase of orbital modulation dependent on the mass loss rate, and subsequently, on the accretion rate after a viscous time of $\sim 10^{5}$ s (Zdziarski et al., 2007b). The variation in mass loss rate could be explained by the triple model as described in Section 1, where the eccentricity variation of the binary system induced by the third companion result in changes to the mass loss rate. Although the superorbital modulation is probably not a consequence of a third companion, the the discovery by Zdziarski et al. (2007b). implies that the superorbital modulation of 4U 1820-30 is due to changes in mass loss rate.

One possible cause of variations in the accretion flow, aside from the presence of a third companion, is irradiation-induced mass transfer instability. This model has ever proposed to explain the flux variation of soft X-ray transients (Hameury et al., 1986) and was included into the hybrid model proposed by Šimon (2003). Due to small binary separation of 4U 1820-30, the irradiation on the companion by the X-ray emission from the neutron star and the inner part of accretion disk is strong. Because the companion of 4U 1820-30 is only partially degenerate (Rappaport et al., 1987), irradiation on the non-degenerate envelope enhances the mass loss of the companion. Chou & Jhang (2023) estimated that at least 40% of the mass lost from the companion is ejected from the binary system. Such a strong outflow is probably caused by the irradiation on the companion, as proposed by Tavani (1991). However, a part of X-ray irradiation on the companion is blocked by the accretion disk, with the area that depending on the scale height of disk rim. When the scale height of disk rim is small, a larger irradiation area enhances the mass loss rate and the accretion flow, which increases the scale height of the accretion disk rim. Conversely, when the scale height of accretion disk rim is large, a larger portion of the companion’s surface is shielded by the disk. This results in a reduction in the mass loss rate, as well as in the accretion flow and the scale height of the accretion disk rim. This a cyclical process may explain the quasi-periodic superorbital modulation of 4U 1820-30.

Suppose the accretion disk in 4U 1820-30 is geometrically thin and optically thick (Pringle, 1981). The shielded region on the companion can be estimated. For 4U 1820-30 with a neutron star mass of $m_{1}=1.4M_{\sun}$, a companion mass of $m_{2}=0.07M_{\sun}$ (Chou & Jhang, 2023), and an orbital period of $P_{orb}=685$ s, we derived the binary separation of $a=1.33\times 10^{10}$ cm from Kepler’s third law. The discovery of the superhump in 4U 1820-30 system (Wang & Chakrabarty, 2010) indicates that the rim of the accretion disk reaches a 3:1 resonance radius, giving a disk radius of $r_{d}=6.4\times 10^{9}$ cm. The scale height of the accretion disk can be evaluated as $\sqrt{kT/\mu m_{H}}/\Omega_{k}$ (Spruit, 2010) where $T$ and $\Omega_{k}$ are the temperature and the Keplerian angular velocity at the disk radius $r$, $m_{H}$ is the mass of a hydrogen atom, $\mu=4$ for a helium-dominated disk, and $k$ is the Boltzmann constant. The temperature at the disk rim is estimated as