© W. Liu et al. 2022

1. Introduction

According to the luminosity class of the optical companion, high-mass X-ray binaries (HMXBs) are divided into supergiant X-ray binaries and Be/X-ray binaries (BeXBs; Reig 2011). Most of the optically identified HMXBs (or HMXB candidates) are known or suspected BeXBs (Liu et al. 2006). The optical companion of a BeXB is a Be star, which is a nonsupergiant, fast-rotating B-type (but it may also include late O-type stars; Negueruela et al. 2004), and luminosity class III-V star that has shown emission lines at some point in its life (Rivinius et al. 2013). There are two different disks in Be/X-ray binaries: circumstellar disks around the Be stars, and accretion disks around the neutron stars, which temporarily appear during X-ray outbursts (Ziolkowski 2002; Hayasaki & Okazaki 2004). BeXBs are classified into persistent and transient sources according to their X-ray properties (Reig & Roche 1999). Transient BeXBs display two types of X-ray outbursts when they are active: type I (or normal) outbursts, and type II (or giant) outbursts. The peak luminosity during type I outbursts is typically L_X ≤ 10³⁷ erg s⁻¹, while during type II outbursts, it may reach the Eddington limit, L_X ∼ 10³⁸ erg s⁻¹.

Swift J0243.6+6124 was detected in X-rays for the first time by the Swift/BAT on 3 October 2017 (Kennea et al. 2017). It is the first Be/X-ray binary emitting at super-Eddington luminosity in our galaxy. Reig et al. (2020) estimated the spectral type and rotational velocity of the companion of Swift J0243.6+6124 to be O9.5Ve and νsini = 210 ± 20 km s⁻¹. X-ray pulsations with a period of ∼9.86 s were detected by Swift/XRT and Fermi/GBM (Jenke & Wilson-Hodge 2017). The orbital period is 28.3 days, and the eccentricity is 0.092 (Doroshenko et al. 2018).

There are several different reported values for the distance to this source in the literature: 2.5 ± 0.5 kpc (Bikmaev et al. 2017) and ∼5 kpc (Reig et al. 2020), both based on optical photometric observations, a lower limit of 5 kpc set by van den Eijnden et al. (2018), ∼6 kpc based on two accretion torque models (Zhang et al. 2019), and $5.5_{-0.3}^{+0.4}$ kpc given in Gaia DR3 catalog (Gaia Collaboration 2016, 2022). If the distance of the source is assumed for 5 kpc, its peak luminosity is estimated as 1 × 10³⁹ erg s⁻¹ (0.1–10 keV). This exceeds the Eddington limit for the neutron star during the giant outburst.

The magnetic field strength of the neutron star in Swift J0243.6+6124 is estimated to be approximately 10¹³ G (Tsygankov et al. 2018; Zhang et al. 2019), although Doroshenko et al. (2018) advocated for a lower value, at the lower limit of the range (3−9)×10¹² G. Based on the detection of electron cyclotron resonance scattering features (CRSFs), Kong et al. (2022) estimated a surface magnetic field of ∼1.6 × 10¹³ G for Swift J0243.6+6124, which unambiguously proves the presence of multipole field components close to the surface of the neutron star. This measured surface magnetic field is the strongest of all known neutron stars with detected electron CRSFs, and it is also the strongest for all neutron star ultraluminous X-ray sources. All types of X-ray binaries have been observed to launch jets, with the exception of neutron stars that have strong magnetic fields (stronger than 10¹² G), which implies that their magnetic field strength restrains jet formation (van den Eijnden et al. 2018). Therefore, the detection of radio emission from Swift J0243.6+6124 during the X-ray outbursts is a surprising result. However, the radio luminosity is two orders of magnitude dimmer than those seen in other accreting neutron stars with similar X-ray luminosities (van den Eijnden et al. 2018), which implies that the magnetic field of neutron stars still plays an important role in the power of launching jets. In this work, we report new optical spectroscopic observations and photometric observations. These observations witnessed the partial dissipation after the giant outburst and subsequent reformation of the Be star’s circumstellar disk.

2. Observations

2.1. Optical spectroscopy

Optical spectroscopic observations were mainly obtained with two telescopes at two different observatories: The observations from the Xinglong Station of National Astronomical Observatories in Hebei province (China) were obtained with the spectrometer OptoMechanics Research (OMR) or BAO Faint Object Spectrograpy and Camera (BFOSC) on the 2.16 m telescope, and the observations from the Lijiang station of Yunnan Astronomical Observatory in Yunnan province (China) used the spectrometer Yunnan Faint Object Spectrograpy and Camera (YFOSC) on the 2.4 m telescope. The OMR was equipped with a 1024 × 1024 (24 micron) pixel TK1024AB2 CCD. The OMR Grism 4 is 1200 lp mm⁻¹, giving a nominal dispersion of 1.2 Å pixel⁻¹, and covering the wavelength ranges 5500–6900 Å. The BFOSC was equipped with a 2048 × 2048 (15 micron) pixel Loral Lick 3 CCD. The dispersion of BFOSC Grism 4 and 8 is 2.97 and 1.20 Å pixel⁻¹, covering the wavelength ranges 4000–8700 and 5800–8280 Å, respectively. The YFOSC was equipped with a 2k × 4k (13.5 micron) pixel E2V 42-90 CCD. The dispersion of YFOSC Grism 8 is 1.47 Å pixel⁻¹, covering the wavelength ranges 4970–9830 Å. In addition, we analyzed new optical spectroscopic observations obtained from the 1.3 m telescope of the Skinakas observatory (SKO) in Crete (Greece). The 1.3 m telescope of the SKO was equipped with a 2048 × 2048 (13.5 micron) pixel ANDOR IKON CCD and a 1302 lines mm⁻¹ grating, giving a nominal dispersion of ∼0.8 Å pixel⁻¹.

We used the Image Reduction and Analysis Facility (IRAF)¹ software package to reduce and analyze all the spectra, performing bias-subtracted correction and flat-field correction on the data, and then removing cosmic rays. A helium-argon calibration lamp was employed to obtain the pixel-wavelength relation. In order to ensure the consistency of spectral processing, all spectra were normalized to adjacent continua. We measured the equivalent width of the Hα lines (hereafter EW(Hα) for short) for five times, each measurement with a different selection of the continuum. The final EW(Hα) is the average of the five measurements, and the error is the standard deviation. The typical error of EW(Hα) is within 5%. The value of the error is determined by the quality of the continuum. The equivalent widths of the He I λ6678 lines (hereafter EW(He I λ6678) for short) were obtained following the same method as for EW(Hα).

The log of the spectroscopic observations is given in Table A.1. EW(He I λ6678) and EW(Hα) are plotted in the third and sixth panels of Fig. 1, respectively. The evolution of the Hα line profiles is plotted in Fig. 2. The evolution of the log(V/R) and the peak separation of Hα line are plotted in Fig. 3 and listed in Table B.1.

Fig. 1. Long-term variations of Swift J0243.6+6124 from 2017 to 2022. First panel: X-ray band light curves. The data from Swift/BAT (15–50 keV) and MAXI (2–20 keV) are labeled with red points and blue downward triangles, respectively. The first day of each year is also marked at the top of the first panel in the format YYYY. Second panel: barycentered and orbit-corrected spin-frequency history measured with Fermi/GBM. Third panel: equivalent widths of He I λ6678 lines. Fourth panel: long-term light curves of the optical V band from AAVSO, ASAS-SN, Reig et al. (2020), and this work. They are separately labeled with light green downward triangles, blue stars, red plus signs, and black points. Fifth panel: long-term light curves of NEOWISE infrared magnitudes at 3.4 μm and 4.6 μm. They are labeled with black plus signs and red points, respectively. Sixth panel: equivalent widths of Hα lines from this work and Reig et al. (2020). They are labeled with black points and red plus signs, respectively. Seventh panel: evolution of the (B − V) color index. The vertical lines separate the three phases into which the observations are divided.

Fig. 2. Evolution of Hα line profiles. The vertical lines mark the rest wavelength of the Hα lines and diffuse interstellar bands (DIB) 6613 Å. All spectra have been normalized with adjacent continua.

Fig. 3. Variations in spectral features of Hα line. Top panel: evolution of log(V/R). Bottom panel: peak separation of Hα line. Stars represent the spectral profiles with V >  R, and triangles represent the spectral profiles with V < R.

2.2. Optical photometry

Optical photometric observations were obtained from five telescopes at three different observatories: From the Xinglong station of the National Astronomical Observatories, Chinese Academy of Sciences (NAOC), observations were obtained with the Tsinghua-NAOC Telescope (TNT, 80 cm), the 60 cm telescope, and the 2.16 m telescope; from the Lijiang station of Yunnan Observatories (YNAO), the data came from the 2.4 m telescope; and from the Yaoan astronomical observation station of Purple Mountain Observatory (PMO), the data came from the Yaoan High Precision Telescope (YAHPT, 80 cm’).

The TNT (80 cm) is an equatorial-mounted Cassegrain system with a focal ratio of f/10, made by AstroOptik, funded by Tsinghua University in 2004 and jointly operated with NAOC, which is equipped with a PI VersArray 1300B LN 1340 × 1300 thin, back-illuminated CCD with a 20 μm pixel⁻¹ size (Wang et al. 2008; Huang et al. 2012). In this configuration, the plate scale is 0.52″pixel⁻¹ and gives a field of view of 11.4 × 11.1 arcmin². The 60 cm telescope is an equatorial-mounted system with a focal ratio of f/4.23, which is equipped with the Andor DU934P-BEX2-DD 1024 × 1024 CCD and provides a field of view of 18 × 18 arcmin². The 2.4 m telescope is an altazimuth-mounted Cassegrain system with a focal ratio of f/8, which is equipped with an E2V CCD42-90 2k × 2k thin, back-illuminated, deep-depletion CCD with a 13.5 μm pixel⁻¹ size. In this configuration, the plate scale is 0.28″pixel⁻¹ and gives a field of view of 9.6 × 9.6 arcmin². The YAHPT (80 cm′) is an altazimuth-mounted, RC optical system with a focal ratio of f/10, made by Astro Systeme Austria, which is equipped with a PIXIS 2048B back-illuminated CCD with a 13.5 μm pixel⁻¹ size. In this configuration, the plate scale is 0.347″pixel⁻¹, providing a field of view of 11.8 × 11.8 arcmin². The 2.16 m telescope is an equatorial-mounted, RC optical system with a focal ratio of f/9, made by NAOC, CAS Nanjing Astronomical Instruments Co., LTD (NAIRC), and the Institute of Automation of the Chinese Academy of Sciences (CASIA), which is equipped with an Andor-DZ936-BEX2-DD 2048 × 2048 CCD with a 13.5 μm pixel⁻¹ size. In this configuration, the plate scale is 0.274″pixel⁻¹, providing a field of view of 9.36 × 9.36 arcmin².

In all five telescopes, Swift J0243.6+6124 was observed through the standard Johnson-Cousins B, V, R, and I filters. The photometric data reduction was performed using standard routines and aperture photometry packages (some from the zphot package) in IRAF, including bias subtraction and flat-field correction. In order to derive the variation in the optical brightness, we selected the reference star Gaia 465628266540345216 (α: 02 43 38.23, and δ: +61 26 40.7, J2000) (according to Reig et al. 2020 the average magnitudes of the reference star are B = 13.67 ± 0.01, V = 13.02 ± 0.01, R = 12.65 ± 0.01, and I = 12.25 ± 0.02) in the field of view of Swift J0243.6+6124 to derive its differential magnitudes. The photometric magnitudes are given in Table C.1.

To study the long-term optical variability of the source, we used the public optical photometric data from the ASAS–SN² Variable Stars Database (Shappee et al. 2014; Jayasinghe et al. 2019). There is a slightly fainter star at 6.2 arcsec from Swift J0243.6+6124 within the full width at half maximum (FWHM). This star is resolved in our photometry. However, the pixel scale and the FWHM in ASAS–SN are 8 arcsec and ∼2 pixels, and hence in these images, the neighboring star contributes to the measured flux from Swift J0243.6+6124. The calibrated V-band magnitude of the fainter close star is V = 14.52 ± 0.01 mag (Reig et al. 2020). We removed the brightness of the neighboring star from the total observed flux. The applied corrections are in the range $\mathrm{\Delta }V=0.{20}_{-0.04}^{+0.03}$ mag. We also made use of the public optical photometric data from the international database of the American Association of Variable Star Observers (AAVSO³). Finally, we also included the optical photometric data from Reig et al. (2020).

The Johnson V-band light curve is plotted in the fourth panel of Fig. 1. A detailed view of the 2017 outburst, including pre- and post-outburst observations, is shown in Fig. 4. The evolution of the (B − V) color index is plotted in the seventh panel of Fig. 1, and the variation of the (B − V) color index versus V-band magnitude is plotted in Fig. 5, where only the data from the 80 cm telescope, the 2.4 m telescope, and the 2.16 m telescope are shown.

Fig. 4. V and IR-band light curves of Swift J0243.6+6124. Top panel: long-term light curves of the optical V band from MJD 56 600 to 59 670. Bottom panel: simultaneous light curves of NEOWISE infrared magnitudes at 3.4 μm and 4.6 μm.

Fig. 5. Variation in (B − V) color index versus V-band magnitude.

2.3. NEOWISE photometry

We made use of the light curves in the W1 (3.4 μm) and W2 (4.6 μm) bands provided by the NEOWISE (Mainzer et al. 2011) project through the IRSA viewer⁴. We plot them in the fifth panel of Fig. 1 and in the bottom panel of Fig. 4.

2.4. X-Ray observations

The Burst Alert Telescope (BAT)⁵ on board Swift (Krimm et al. 2013), MAXI⁶, and the Gamma-ray Burst Monitor (GBM)⁷ on board Fermi (Meegan et al. 2009) have been monitoring Swift J0243.6+6124 in the hard X-ray energy band (15–50 keV with BAT, 2–20 keV with MAXI, and 12–50 keV with GBM) since October 2017. One type II X-ray outburst and several type I outbursts were detected between October 2017 and January 2019. The X-ray band light curves from BAT (15–50 keV) and MAXI (2–20 keV) are plotted in the first panel of Fig. 1. The spin-frequency history measured with GBM is plotted in the second panel of Fig. 1.

3. Results

Figure 1 shows the X-ray, optical, and IR long-term variability of Swift J0243.6+6124. The observations cover the October 2017 X-ray outburst and the changes experienced by the source during the following four-year period. After the X-ray outbursts, the optical brightness of the source decreased and reached a minimum on ∼ MJD 58530. Since then, it has been recovering. There is a clear correlation between the optical and IR flux on long timescales (weeks or months). When the source is active in the X-ray, the optical and IR also correlate with the overall X-ray flux in the sense that the source is bright in the optical and IR at the time of the outbursts. The strength of the Hα line also follows the same general trend, although the minimum after the X-ray outbursts appears to be delayed by ∼100–200 days with respect to the optical or IR continuum flux.

3.1. Three phases of our data

After examining the long-term light curves and spectral evolution, we divided the observations into three different epochs or phases (Fig. 1). These phases reflect significant changes in the properties of the data. Each phase is characterized by a different pattern of X-ray, optical, or IR variability. Phase I corresponds to the giant 2017 X-ray outburst; during phase II, the source experiences a gradual fading of its brightness and a weakening of the spectral line parameters; and in phase III, the long-term trend is reversed and the source exhibits a gradual increase in the brightness and strength of the spectral lines, most significantly, in the Hα line.

Phase I (MJD 58030–58180) covers the giant X-ray outburst. The X-ray luminosity changes by about two orders of magnitude (10³⁷ − 10³⁹ erg s⁻¹; Doroshenko et al. 2020) and the V band by 0.4 mag. EW(Hα) and EW(He I λ6678) display erratic variability. The (B − V) color displays the largest variation in the entire period of the observations, with a change of about 0.1 mag in 25 days.

In phase II (MJD 58180–58530), the X-ray variability is characterized by regular type I outbursts with changes in luminosity of about one order of magnitude between 10³⁶ − 10³⁷ erg s⁻¹. The brightness in the V band decreases by 0.1 mag and in the near-infrared by 0.6 mag. EW(Hα) decreases from a maximum of −11 Å to a minimum of −5 Å. EW(He I λ6678) presents large scatter, but it also decreases on average. The dispersion in EW(He I λ6678) measurements is most likely due to the low signal-to-noise ratio (S/N). There are no suitable B-band observations during the first half of this phase; thus we do not have (B − V) data. During the second half, the (B − V) color is lower on average (i.e., bluer emission) than during phase I.

In phase III (MJD 58530–), the source is no longer detected in X-ray, while all optical and infrared indicators increase gradually. At the end of this phase, EW(Hα) and the V-band magnitude recover to almost pre-outburst values. The overall optical emission becomes redder as the system evolves in this phase.

3.2. Hα line profiles

Figure 2 shows the profiles of the Hα line. Although we are limited by the low spectral resolution, some general variability trends are visible. The Hα emission lines present asymmetrically blue-dominated profiles (V >  R) between October 2017 and September 2018 (phase I–phase II). The red peak is only noticeable in the earliest spectrum on 27 October 2017. The central rest wavelength (marked by the vertical dashed line in Fig. 2) lies systematically to the right of the peak, confirming the asymmetry. At the end of 2018 and during 2019 (phase II–phase III), the line flux gradually shifts toward the red. Approximately symmetric single-peak profiles are observed in the 2020–2021 spectra (phase III). At the end of 2021, the profiles are again blue-dominated.

The He I line also exhibits a variable asymmetric double-peak profile in phase and with intensity variations similar to those of the Hα line. Because the He I line forms much closer to the Be star than the Hα line, the He I V/R variability implies that the density perturbation affects the inner parts of the disk as well as the outer parts where the Hα line is formed. Unfortunately, the low S/N and low spectral resolution prevented a detailed analysis of this line.

4. Discussion

We have been monitoring the optical counterpart of the Be/X-ray binary Swift J0243.6+6124 since its discovery in October 2017 and obtained optical spectra on a regular basis. Using data from both the archives and the literature, we are able to study the long-term variability over at least four years after the outburst. We attribute this long-term variability to the evolution of the Be star’s circumstellar disk.

The variability patterns of the X-ray, optical, and IR parameters allowed us to divide the observations into three different phases. Phase I covers the rise and decay of the X-ray outburst. The largest and fastest variability timescales are seen during this phase. This violent event must have affected the entire structure of the circumstellar disk. In particular, the Hα line displayed several sudden changes of 10%–20% in its strength on timescales of about one month during this phase (see the sixth panel of Fig. 1). These fast changes could be attributed to an inhomogeneous warped disk (Reig et al. 2020). Alternatively, because the fast changes of EW(Hα) appear to be modulated with the orbital period of the system (P_orb = 28.3 days), reprocessed emission might be at their origin. The high spin-up rate (as seen in the second panel of Fig. 1) suggests that a short-lived accretion disk also formed around the neutron star.

Phase II corresponds to a low optical state in which not only the optical and IR continuum flux gradually decreased, but the Hα emission line also became weaker. Type II outbursts involve the accretion of a large amount of material from the equatorial disk and almost always lead to the dissipation of the entire or part of the disk (Reig et al. 2016). We identify this phase with the dissipation of the disk. The evolution of the X-ray and optical emissions provides some clues about how this dissipation took place. The presence of type I outbursts indicates that the neutron star accretes material from the disk in subsequent periastron passages, while at the same time, the decrease in the (B − V) color indicates that the emission becomes bluer, as expected from a smaller or more compact disk. These results can be understood if the outer parts of the disk cease to be bound, that is, are expelled, while the inner parts of the disk collapse toward the star once the disk formation mechanism stops. The higher X-ray intensity of the last type I outburst in this phase may be due to the fact that the neutron star encountered a higher-density part of this distorted disk. This phase ends when the magnitude and colors reach a minimum. The system did not lose the disk entirely, as the Hα line remained with an emission profile and the minimum EW(Hα) was still ∼ − 5 Å. Even the He I λ6678 did not revert into an absorption profile. The equivalent width of this line remained ∼0 Å during 2020 and 2021, indicating that the line was filled with emission.

Phase III represents a phase in which the disk grew again. The optical and IR flux increased, as did the strength of the Hα line. The overall emission became redder (i.e., (B − V) increased). The latest observations show that EW(Hα) and EW(He I λ6678) approached pre-outburst values. After an initial brightening period that lasted for about a year, the V-band magnitude stabilized at a level of 12.8^m, at which it has remained since September 2020.

Photometric variability with a characteristic period of ∼1250 days (MJD 40000–41250) and an amplitude of ∼0.15 mag in the B-band was reported by Nesci (2017) based on archival data from the Asiago Observatory taken between 1967 and 1976. This timescale is similar to the one during MJD 57250 to 58530 (see Fig. 4). When the V-band flare is ignored, the underlying trend of the V-band observations correlates very well with the IR-band observations. Given the strength of this correlation, the lack of a similar flare in the IR light curve can be attributed to the low cadence of the IR observations. The disk starts to grow on ∼ MJD 57250, reaching a maximum size on ∼ MJD 58000 when it begins to decline, and reaching a minimum on ∼ MJD 58530. The variation in the V/IR-band can be interpreted in terms of the evolution of the Be star’s disk. This long-term smooth change in brightness is due to the formation or growth of the disk and its subsequent dissipation. We estimate that the overall timescale in phase III for the formation and dissipation of the circumstellar disk is about 1500 days. The dissipation phase is significantly faster, about 300 days, than the formation phase of about 1200 days.

Although the evolution of the optical and IR parameters is affected by observational gaps, because the position of the source is too close to the Sun, Fig. 1 shows a delay between the minima of the optical and IR continuum (which marks the shift from phase II to phase III) and the minima of EW(Hα). This delay may be understood by invoking the different sites of the continuum and discrete emission in the equatorial disk of a Be star. According to Carciofi (2011), the disk V band is typically formed very close to the star, within about 2 R_star. In contrast, the Hα emission line is formed at larger radii (Slettebak et al. 1992). If this interpretation is correct, then the minimum flux detected in the continuum first would imply that the dissipation of the disk began from the inner parts. The type I regular X-ray outbursts during phase II indicate that the accretion of the neutron star persisted for about one year after the main outburst. In terms of the strength of the Hα line, the giant X-ray outburst took place when EW(Hα) ∼ − 11 Å. The increase rate of EW(Hα) during phase III implies that the source will reach this value again on ∼MJD 60250. If we take this value (−11 Å) as the triggering value of the giant outburst, then we should expect another large event by the end of 2023.

The He I λ6678 follows the same long-term trend as the Hα line: fast and large amplitude changes during phase I, a weakening during phase II, and a slow recovery during phase III. The latest spectrum in October 2021, in which a small peak started to develop and EW(He I λ6678) < −0.1 Å, marks the formation of a new emission line. We note that the giant X-ray outburst occurred when EW(He I λ6678) ∼ −(0.2−0.5) Å.

4.1. Hα line profile variability and V/R ratio

The V/R variability is defined as the intensity variations of the two peaks (known as violet and red peaks) in the split profile of a spectral line. In many Be stars, if they are monitored over a long enough period of time, these variations are quasi-periodic (Okazaki 1997). We define the V/R ratio of the Hα line as V/R = (I(V)−I_c)/(I(R)−I_c), where I(V), I(R), and I_c are the intensities of the violet peak, red peak, and continuum, respectively. We also measured the separation of the violet and red peaks by fitting two Gaussian functions to the spectral line profile. When the disk velocity is assumed to be Keplerian, the peak separation gives a measure of the velocity field. There is no obvious trend in the peak separation between different spectral line profiles, which is mainly distributed around 175–250 km s⁻¹. The V/R ratios and the peak separation of the Hα line are listed in Table B.1 and plotted in Fig. 3.

The V/R variability has been associated with density perturbations in the disk (Hanuschik et al. 1995). When this density perturbation moves around inside the disk, the profile changes. We observe a blue-dominated profile (V >  R) in 2017 that turned into an almost single peak profile (V ∼ R) in 2018 (see Fig. 4 of Reig et al. 2020), and a red-dominated profile (V <  R) from the end of 2018 to 2019. The spectra in 2021 return to blue-dominated profiles (V >  R). Thus, we may have covered an entire V/R cycle. The V/R quasi-period would be about four years, which is normal for BeXBs (Mennickent et al. 1997).

In principle, the question of whether the motion of the perturbation occurs in the same sense (prograde rotation) or opposite sense (retrograde rotation) to the stellar rotation can be determined from the observations. Telting et al. (1994) realized that a prograde rotation implies that (I) a V >  R phase, (II) a shell absorption profile, (III) a V <  R phase, and (IV) a weak central absorption profile will appear in order. A retrograde rotation would give rise to the reversed sequence: (IV)→(III)→(II)→(I). Because of the small disk inclination (Reig et al. 2020), we cannot distinguish a prograde or retrograde rotation in the characteristic line shapes. These characteristic line shapes can translate into noticeable photometric variations, however. According to Mennickent et al. (1997), we can expect a minimum brightness when V = R prior to the V <  R (V >  R) phase if the motion is prograde (retrograde). In Swift J0243.6+6124, the minimum brightness in the photometric V band occurred during the V = R phase before the V <  R phase began, ∼MJD 58450, confirming the prograde nature of the precession inside the disk.

4.2. Variation in (B − V) color index and inclination of the Be star’s disk

Figure 5 shows the (B − V) color index as a function of the V magnitude. We mark the different variability phases defined above with different colors. It has been noted that this kind of plot can be used to constrain the inclination angle of the system (Harmanec 1983). Systems that show a positive correlation, that is, as the disk forms (or equivalently, as EW(Hα) increases), the optical intensity increases and the emission becomes redder (i.e., (B − V) increases), are thought to be seen at small or moderate inclination angles, while systems that show a negative correlation, in which the optical intensity decreases even though the disk is growing (EW(Hα) and (B − V) increase), are associated with large inclination angles.

Harmanec (1983) introduced the concept of a pseudophotosphere to explain this effect. At large inclination angles (for equator-on stars), the inner parts of the Be envelope partly block the stellar photosphere, and thus the optical brightness decreases. Meanwhile, the overall emission becomes redder because the contribution of the disk increases. At small or intermediate inclination, as the disk grows, an overall (star plus disk) increase in brightness is expected.

Figure 5 suggests that Swift J0243.6+6124 is viewed at small or intermediate angles. This result is consistent with the 30° angle estimated by Reig et al. (2020) from the emission-line profile.

5. Conclusions

We have conducted spectroscopic and photometric observations of Swift J0243.6+6124 to study the changes in the structure of the circumstellar disk of the Be star and the material transfer between the Be star’s circumstellar disk and the neutron star during 2017–2022. This period covers a giant X-ray outburst (type II) and several orbit-modulated outbursts (type I). We divided our data into three phases based on the intensity of the X-ray, optical, and infrared emission. Phase I covers the 2017 X-ray outburst. In this phase, the source reaches the largest equivalent width of the Hα and He I λ6678 lines and the brightest X-ray, optical, and infrared intensity, characterized by outbursts or flares. In phase II, the source displays a long-term decrease in the optical and infrared intensities. The X-ray band is characterized by the occurrence of several minor orbit-modulated outbursts. In phase III, the source exhibits a long-term brightening of the optical and infrared magnitudes, and the equivalent width of the Hα and He I λ6678 lines also increases. During this phase, no X-ray activity is observed. We focused on the correlated optical, IR, and X-ray long-term variability. During the period of our observations following the giant X-ray outburst, the optical and IR continuum flux and the strength of the Hα line first decreased and then increased. We interpret this long-term variability in terms of the dissipation and reformation of the Be star’s circumstellar disk. Although the 2017 giant X-ray outburst is the most luminous ever recorded in a BeXB, it does not lead to the complete loss of the Be star’s disk. We estimate that at the rate at which the disk is reforming, the equivalent width of the Hα line will reach pre-outburst values in about 1–1.5 years.

Acknowledgments

We acknowledge the support of the staff of the Xinglong 2.16 m telescope, the Xinglong 80 cm telescope, and the Xinglong 60 cm telescope. This work was partially supported by the Open Project Program of the CAS Key Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese Academy of Sciences. We acknowledge the support of the staff of the Lijiang 2.4 m telescope. Funding for the telescope has been provided by CAS and the People’s Government of Yunnan Province. Skinakas Observatory is run by the University of Crete and the Foundation for Research and Technology-Hellas. This research has made use of data provided by the Yaoan High Precision Telescope. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. Swift/BAT transient monitor results provided by the Swift/BAT team. Fermi/GBM results provided by the Fermi Science Support Center. This publication makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration. This work is supported by the National Key R&D Program of China (2021YFA0718500) and the National Natural Science Foundation of China (Grants No. U2031205, 11733009). X. Wang is supported by the National Science Foundation of China (NSFC grants 12033003 and 11633002), the Scholar Program of Beijing Academy of Science and Technology (DZ:BS202002), and the Tencent Xplorer Prize.

References

Appendix A: Table of spectroscopic observations

Appendix B: Table of V/R ratio and peak separation

Appendix C: Table of photometric observations

All Tables

All Figures

Fig. 1. Long-term variations of Swift J0243.6+6124 from 2017 to 2022. First panel: X-ray band light curves. The data from Swift/BAT (15–50 keV) and MAXI (2–20 keV) are labeled with red points and blue downward triangles, respectively. The first day of each year is also marked at the top of the first panel in the format YYYY. Second panel: barycentered and orbit-corrected spin-frequency history measured with Fermi/GBM. Third panel: equivalent widths of He I λ6678 lines. Fourth panel: long-term light curves of the optical V band from AAVSO, ASAS-SN, Reig et al. (2020), and this work. They are separately labeled with light green downward triangles, blue stars, red plus signs, and black points. Fifth panel: long-term light curves of NEOWISE infrared magnitudes at 3.4 μm and 4.6 μm. They are labeled with black plus signs and red points, respectively. Sixth panel: equivalent widths of Hα lines from this work and Reig et al. (2020). They are labeled with black points and red plus signs, respectively. Seventh panel: evolution of the (B − V) color index. The vertical lines separate the three phases into which the observations are divided.

In the text

	Fig. 2. Evolution of Hα line profiles. The vertical lines mark the rest wavelength of the Hα lines and diffuse interstellar bands (DIB) 6613 Å. All spectra have been normalized with adjacent continua.
In the text

	Fig. 3. Variations in spectral features of Hα line. Top panel: evolution of log(V/R). Bottom panel: peak separation of Hα line. Stars represent the spectral profiles with V >  R, and triangles represent the spectral profiles with V < R.
In the text

	Fig. 4. V and IR-band light curves of Swift J0243.6+6124. Top panel: long-term light curves of the optical V band from MJD 56 600 to 59 670. Bottom panel: simultaneous light curves of NEOWISE infrared magnitudes at 3.4 μm and 4.6 μm.
In the text

	Fig. 5. Variation in (B − V) color index versus V-band magnitude.
In the text