Long-term Study of the First Galactic Ultraluminous X-Ray Source Swift J0243.6+6124 Using NICER

We present the results obtained from detailed X-ray timing and spectral studies of X-ray pulsar Swift J0243.6+6124 during its giant and normal X-ray outbursts between 2017 and 2023 observed by the Neutron star Interior Composition Explorer (NICER). We focused on a timing analysis of the normal outbursts. A distinct break is found in the power density spectra of the source. The corresponding break frequency and slopes of the power laws around the break vary with luminosity, indicating a change in the accretion dynamics with the mass accretion rate. Interestingly, we detected quasiperiodic oscillations within a specific luminosity range, providing further insights into the underlying physical processes. We also studied the neutron star spin period evolution and a luminosity variation in the pulse profile during the recent 2023 outburst. The spectral analysis was conducted comprehensively for the giant and all other normal outbursts. We identified a double transition at luminosities of ≈7.5 × 1037 and 2.1 × 1038 erg s−1 in the evolution of continuum parameters like the photon index and cutoff energy with luminosity. This indicates three distinct accretion modes experienced by the source, mainly during the giant X-ray outburst. A soft blackbody component with a temperature of 0.08–0.7 keV is also detected in the spectra. The observed temperature undergoes a discontinuous transition when the pulsar evolves from a sub- to super-Eddington state. Notably, in addition to an evolving 6–7 keV iron line complex, a 1 keV emission line was observed during the super-Eddington state of the source, implying X-ray reflection from the accretion disk or outflow material.


INTRODUCTION
Ultraluminous X-ray sources (ULXs) are non-nuclear point-like X-ray sources observed in extra-galactic regions that exhibit luminosities exceeding 10 39 erg s −1 (for a detailed review see Kaaret et al. 2017).Initially, these X-ray sources were speculated to be intermediatemass black holes (IMBH) (Colbert & Mushotzky 1999).Later, it was observed that the ULXs are stellar-mass black holes or neutron stars that are present in close binaries in which the donor fills the Roche lobe and accretion occurs in supercritical/super-Eddington mode (King et al. 2001).In 2014, the first Ultraluminous Xray pulsar (ULXP) was discovered, with coherent pulsations detected from the ULX M82 X-2 by Bachetti et al. (2014).This finding provided a unique opportunity to study super-Eddington accretion onto the magnetized neutron stars.Since then, several ULXPs have been detected, contributing to our understanding of these systems (Mushtukov et al. 2019 & references therein).
The accretion process and radiation conversion in ULXPs and typical accretion-powered X-ray pulsars share similarities.Accretion-powered X-ray pulsars possess a strong magnetic field that directs the infalling matter toward the polar regions of the neutron star.Due to the intense gravitational field near the neutron star, the accreting material gains substantial kinetic energy.Eventually, depending on the mass accretion rate, the material decelerates through the gas or radiative shocks as it reaches the neutron star's surface.The combination of accumulated material and generated radiation in that region forms an accretion column, which serves as the dominant source of radiation in accreting pulsars (Becker & Wolff 2007;Becker et al. 2012;Reig & Nespoli 2013).The primary emission in accreting pulsars consists of seed photons, in soft X-rays, originating from the hot spots above the magnetic poles.These photons undergo Compton upscattering with the infalling plasma giving rise to broadband X-ray emis-sion (Becker & Wolff 2007).The broadband continuum spectra are typically described by a blackbody component and a cutoff power-law component (Becker & Wolff 2007;Reig & Nespoli 2013).These systems also exhibit emission features mostly from iron (Fe) and absorption features in the 10-100 keV range caused by the cyclotron resonant scattering (Walter et al. 2015;Jaisawal & Naik 2016;Staubert et al. 2019;Chhotaray et al. 2023).The latter provides a direct means to estimate the magnetic field strength of the pulsar.
In addition to their spectral characteristics, accreting pulsars display various temporal properties.They exhibit coherent pulsations corresponding to their spin periods and often undergo spin-up phases, where the spin-up rate increases with increasing luminosity (Lamb et al. 1973).The pulse profiles of accreting pulsars also show variation with the change in mass accretion rate onto the neutron star (Nagase 1989;Wilson-Hodge et al. 2018;Jaisawal et al. 2023).The power density spectrum (PDS) of pulsars exhibits intriguing characteristics.It typically exhibits narrow features, indicating the presence of pulsations, as well as broad features that are attributed to quasi-periodic oscillations (QPOs) in the source (van der Klis et al. 1987;Raichur & Paul 2008;Wilson-Hodge et al. 2018).
A remarkable transient X-ray outburst observed between 2017 and 2018 led to the discovery of Swift J0243.6+6124.Historically, the Swift/BAT (15-50 keV) telescope detected this source at a flux level of approximately 80 mCrab on October 3, 2017, marking the onset of its 2017-2018 giant X-ray outburst (Kennea et al. 2017).Swift J0243.6+6124 is considered to be the first Galactic ULX due to its intense X-ray luminosity reaching up to an order of 10 39 erg s −1 (Wilson-Hodge et al. 2018;Jaisawal et al. 2019;Doroshenko et al. 2020).Timing investigations revealed that Swift J0243.6+6124hosts a neutron star with a pulsation period of 9.8 seconds (Kennea et al. 2017;Jenke & Wilson-Hodge 2017;Jaisawal et al. 2018).Optical spectroscopic observations identified the companion (donor) star as an O9.5Ve type star (Kouroubatzakis et al. 2017;Reig et al. 2020).The system is known to have a relatively short orbital period (P orb ) of around 28 days and a mildly eccentric orbit with an eccentricity (e) of approximately 0.1 (Doroshenko et al. 2018;Wilson-Hodge et al. 2018).
A rapid spin-up rate of the pulsar was observed during the giant outburst (Doroshenko et al. 2018).The pulse profile and the pulsed fraction (P F ) exhibited complex variation with luminosity and showed a significant change around the critical luminosity of ∼ 10 38 erg s −1 (Wilson-Hodge et al. 2018;Tsygankov et al. 2018).
A QPO-like feature at 50-70 mHz is reported only in a particular luminosity range in the power density spectra (PDS) (Wilson-Hodge et al. 2018).Further timing analysis by Doroshenko et al. (2020) revealed major changes in pulse profiles and power spectrum at specific luminosity levels.
The broadband continuum of the source during the outburst is well described by an absorbed cutoff power law and black body components (Jaisawal et al. 2018).A two-component transition in the values of the spectral parameters is observed during the giant outburst (Kong et al. 2020).Recent studies provided significant advancements in understanding the magnetic field of this pulsar.Jaisawal et al. (2019) conducted a detailed investigation of the iron line width evolution with luminosity, suggesting a possible origin from the disc would require a magnetic field range of 10 11 to 10 12 Gauss.Additionally, no evidence of a cyclotron absorption line below 100 keV was found (Jaisawal et al. 2018;Beri et al. 2021).Kong et al. (2022) reported the discovery of the CRSF using data from Insight-HXMT in the 120-140 keV range, estimated the magnetic field of the pulsar to be ∼1.6×10 13 Gauss.
Swift J0243.6+6124showed giant and subsequent normal outbursts within MJD 58029 (2017-10-03) andMJD 58533 (2019-02-19).The source also recently went into an outburst phase between MJD 60097 (2023-06-02) and 60190 (2023-09-03) (Figure 1).In this paper, we used Neutron star Interior Composition Explorer (NICER) observations for the long-term study of Swift J0243.6+6124 in the soft X-ray band.The detailed spectral characteristics of the source in the soft X-ray band are still a matter of investigation.Thus, a thorough spectral study of the source during the pulsar's giant and subsequent normal X-ray outbursts, including the new outburst in 2023 is performed in this paper.Moreover, our timing analysis is primarily focused on the multiple normal outbursts, including the recent outburst in 2023, during the post-giant outburst period between MJD 58303-60190.The long-term monitoring capability of NICER allows us to understand the source's spectral and timing characteristics over many magnitudes of luminosity.This paper is organized as follows: Section 2 provides an overview of the observations and data reduction procedures applied to the NICER data.Section 3 presents the results of the timing analysis.Section 4 focuses on results obtained from the spectral analysis.The discussions and conclusions are presented in Sections 5 and 6, respectively.NICER, launched in June 2017 and installed on the International Space Station, is equipped with the X-ray Timing Instrument (XTI; Gendreau et al. ( 2016)) designed to operate in the 0.2-12 keV energy range.The XTI comprises 56 X-ray concentrator optics, each paired with a silicon drift detector, providing non-imaging observations (Prigozhin et al. 2016).These concentrator optics are comprised of 24 nested grazing-incidence goldcoated aluminum foil mirrors of parabolic shape.The XTI offers a high time resolution of approximately 100 ns (rms) and a spectral resolution of about 85 eV at 1 keV.Its field of view covers approximately 30 arcmin 2 in the sky.The effective area of NICER is approximately 1900 cm 2 at 1.5 keV, utilizing 52 active detectors.The XTI is divided into seven groups of eight focal plane modules (FPMs), with each group managed by a Measurement/Power Unit (MPU) slice.This arrangement enables independent operation and control of the FPMs within each group.
We used publicly available NICER data of Swift J0243.6+6124observed between MJD 58029 (2017-10-03) and 58533 (2019-02-19) in our study for giant and subsequent normal X-ray outbursts.These data are stored under observation ids 1050390xxx with a net exposure time of 408 ks.Moreover, for the recent 2023 outburst, we accumulated around 104 ks net exposure data observed from June to September 2023 under observation ids 6050390227-6050390277.We used the nicerl2 pipeline available under HEASoft version 6.30 to process unfiltered event data of NICER.The analysis is performed in the presence of gain and calibration database files of version xti20221001.Good time intervals (GTI) are selected by using nimaketime tool alongside the application of standard filtering criteria based on the South Atlantic Anomaly (SAA), the elevation angle of > 15 • from the Earth limb, a > 30 • offset from the bright Earth, and pointing offset of 0.015 • , on the clean data after nicerl2.The magnetic cutoff rigidity of 4 GeV/c is also applied to filter any possible high energy charge particle backgrounds.We also considered SUNSHINE==0 filtering in GTI to extract the night side data from the 2023 outburst (Ng et al. 2023).This is included due to the known visible-light leak1 in the XTI optical bench of NICER on May 22, 2023.To account for the effects of Earth and satellite motion during the observations, barycentric correction is applied using solar system ephemeris JPL-DE430.We extracted the light curves and spectrum from each observation in XSELECT.The dead time correction to the light curves is not applied as the count rate is below 20000 counts s −12 during these normal outbursts.However, we considered the effect of dead time on spectra during the 2017-18 giant outburst by following Wilson-Hodge et al. (2018).The corresponding spectral background of each observation is obtained using the nibackgen3C50 tool (Remillard et al. 2022).The spectral response matrix and ancillary response files are created using nicerrmf and nicerarf commands, respectively.
Furthermore, in conjunction with NICER, we incorporated the daily X-ray monitoring data obtained from Swift/BAT.These observations encompassed the energy range of 15-50 keV, as detailed in the study conducted by Krimm et al. (2013).To ensure data integrity, the rows marked as good quality data, as indicated by a DATA FLAG value of 0, were included in our analysis.

TIMING ANALYSIS AND RESULTS
We conducted timing analysis during the post-giant outbursts including the 2023 outburst with NICER.We selected these post-giant periods (MJD 58303-60190) to examine the periodic and quasi-periodic oscillations from the neutron star.The timing properties of the pulsar during the 2017-2018 giant outburst are reported in Wilson-Hodge et al. (2018) using NICER data.Figure 1 illustrates the Swift/BAT monitoring light curve (solid red circles) in the 15-50 keV band, and shaded parts above it represent probed regions using NICER data.
The NICER light curves were generated using a bin size of 0.1 seconds in the 0.5-10.0keV energy range.The timing analysis is performed on a total of 105 Ids, each with a minimum good exposure time of 1400s.We searched for the pulsating signal in the light curve following the χ 2 -maximization technique (Leahy 1987) using the efsearch task of FTOOLS package.We also applied the orbital correction to data using the orbital parameters provided by Fermi/GBM team 3 .The orbital correction was done to obtain the intrinsic spin period of the neutron star that gets affected by the binary orbital modulation.In this paper, we present the spin frequency evolution of the source during the recent 2023 outburst (Figure 2).The spin frequency evolution during previous outbursts such as giant and subsequent normal outbursts between MJD 58029 and 58533 is reported in Wilson-Hodge et al. (2018) and Serim et al. (2023), respectively.
The bottom panel of Figure 2 shows the pulse frequency of the pulsar obtained from NICER and publicly 2 https://heasarc.gsfc.nasa.gov/docs/nicer/dataanalysis/nicer analysis tips.h  between MJD 58010 (2017-09-17) and60220 (2023-10-03).The light curve is binned by 1 day time frame.The outbursts that are studied using NICER observations are represented with shaded regions.The dark-violet, violet, and green color mark giant, subsequent normal outbursts, and the 2023 normal outburst phase of the source, respectively.available Fermi/GBM data from the 2023 outburst.We have presented the evolution of 12-50 keV pulsed flux from Fermi/GBM as well as the source luminosity in 0.7-10 keV observed by NICER in the top panel of the figure.The spin frequency during the recent outburst changed from ≈ 102.03 to 102.12 mHz in both NICER and Fermi/GBM data.
We also generated the pulse profile of the pulsar during its normal outbursts including the 2023 outburst to study the geometry of pulsed beamed emission in soft X-rays.For this, each light curve is folded at its corresponding spin period using the efold task of FTOOLS package.The pulse profile variation with luminosity from previous giant and subsequent normal outbursts observed between MJD 58029 and 58533 studied by Wilson-Hodge et al. (2018), Serim et al. (2023), andLiu et al. (2023).We noticed no significant difference in pulse profile evolution with luminosity during the normal outbursts between MJD 58303-58533 and 2023 normal outbursts.Hence, in Figure 3, we show pulse profile evolution with luminosity during the normal outbursts between MJD 58303-58533 and 2023 normal outbursts between MJD 58303-60190 for completeness.The pulse profiles appear complicated with multiple dips/notches at certain pulse phases at luminosities below 6×10 37 erg s −1 .To show the detailed variation of dips/notches and profile evolution with luminosity, we represented the pulse profiles from 12 NICER observations in Fig- ure 4. Below the luminosity of ∼0.5×10 37 erg s −1 , the pulse profiles are single peak dominated in nature.Various dip-like structures arise in the pulse phase range of 0.5-1.0 for luminosity between ∼(0.5-2.5)×10 37erg s −1 .One of the dips became deeper, and the profiles evolved into a double-peaked structure at a luminosity ∼2.5×10 37 erg s −1 .Further, a smooth single peaked profile is observed at luminosities beyond ∼6×10 37 erg s −1 .
Furthermore, we computed the pulsed fraction (P F ) of the pulse profiles from the 2023 outburst to under-  stand the evolution of soft X-ray emission of the pulsar.Such a study has not been performed using NICER, except for the 2017-2018 giant outburst.Thus, we also studied the remaining post-giant outburst data in our study.We calculated the P F using the root mean squared method as given by Wilson-Hodge et al. (2018) (see also Ferrigno et al. 2023).It is expressed as: The variation of P F with luminosity is presented in the panel (a) of Figure 5.We fitted a curve to the obtained P F values at different luminosities using the spline interpolation method, to understand the trend of P F evolution.The red line is the best-fit spline curve, and the shaded area indicates the moving average standard deviation of data points.The P F varied between ∼10% to 15% for the source luminosity in the range of 2 × 10 36 to 9 × 10 37 erg s −1 .A relatively higher value of P F (∼20-25%) was detected at a luminosity of ∼2 × 10 37 erg s −1 .

Power Density Spectrum (PDS) Analysis
The evolution of the PDS and its features during normal outbursts including the 2023 outburst is investigated for the first time in this work.We performed PDS analysis on 105 light curves in the 0.5-10 keV energy range using NICER data to follow the luminosity dependency variation of the PDS.Additionally, we searched for the sign of any QPO in the PDS to compare these features with the same observed during the giant outburst (Wilson-Hodge et al. 2018).
We used the powspec tool from the XRONOS package to generate the power density spectra.The light curves were divided into intervals of ∼400 s, and then the PDS of each interval was generated.The final PDS was obtained by averaging the PDS of segmented light curves.This method improves the detection probability of QPO-like features.We also employed the powspec norm=-2 command to obtain white-noise subtracted averaged PDS.Following this, the power is expressed in units of (RMS/mean) 2 /Hz.The resulting PDS for the observation ID 1050390170 (MJD 58322) is illustrated in Figure 6.
The PDS exhibited narrow peaks at multiple of a main frequency around 0.101 Hz.These peaks correspond to the pulsar's spin frequency and its harmonics, which are ignored during the fitting of PDS (James et al. 2010).Initially, we attempted to fit the PDS continuum by a simple power law model (powerlaw) using XSPEC.However, this model proved inadequate in fitting the overall PDS across a wide frequency range, spanning from approximately 0.001 Hz to 5.0 Hz.Subsequently, we replaced the power law with the broken power-law model (bknpower), which resulted in an acceptable chi-square (χ 2 ) value.We studied the evolution of the slope of power-law before break (Γ 1 ), break frequency (Br f ), and slope of power-law after break (Γ 2 ) in the PDS w.r.t luminosity, which is presented in the panel (d), (e), and (f) of Figure 5, respectively.The red line represents the best-fitted curve using the spline interpolation method to examine the evolution of the parameter.From the  figure, the break frequency Br f varied between 50 and 80 mHz below 7.5×10 37 erg s −1 .Beyond this luminosity, the Br f increases to 140 mHz, surpassing the pulsar frequency (101 mHz).Γ 1 and Γ 2 values varied between 0.5-1 and 2-3, respectively, within the studied luminosity range.
After fitting the PDS continuum with a broken powerlaw model, we examined the residuals for the presence of potential QPOs.In addition to the neutron star's spin frequency and its harmonic components, we observed a broad hump-like residual below the pulsar spin frequency in a few cases (around 16 IDs).To explore this further, we employed a Lorentzian function (lorentz) as the shape of this hump appeared asymmetric.The Lorentzian function has been widely employed in several QPO studies (see e.g., Belloni et al. 2002).Furthermore, the significance of the QPOs is determined using the method described in Boutelier et al. (2010) based on the Lorentzian fitting.We found that QPOs from 10 out of 16 IDs between MJD 58322-58521 exhibited a significance of more than 3σ using the above method.The frequency (QPO f ) and full width at half maximum (QPO W ) of these detected QPOs and their respective significance are presented in Table 1.The evolution of QPO f and QPO W with luminosity is also presented in panels (b) and (c) of Figure 5, respectively.We did not observe any QPO-like feature during the recent 2023 outburst of Swift J0243.6+6124.

SPECTRAL ANALYSIS AND RESULTS
We conducted the spectral analysis using NICER data across multiple outbursts observed between 2017 and 2023.This allowed us to probe changes in the spectral parameters over a broad range of luminosities.The spectral analysis was performed in the energy range of 0.7-10.0keV using the XSPEC (v-12.11.0,Arnaud 1996) package.The above energy range is selected to avoid the spectral uncertainties below 0.4 keV and above 10.0 keV, and a strong edge-like feature appears near 0.5 keV, especially in the brighter observations.We found the 0.5 keV feature depends on the choice of the photo- electric absorption model as well as assumed abundances up to some extent.The 0.5 keV feature may also have a calibration origin from Oxygen edge4 .A systematic uncertainty of 1.5% was also applied, as recommended by the instrument team.For quantifying line-of-sight X-ray absorption, wilm abundance table (Wilms et al. 2000) is used with Vern5 photo-ionization cross section.
The spectra were binned with a minimum of 30 counts per energy bin to allow for the application of chi-square statistics in the analysis.
To understand the evolution of the pulsar emission over multiple outbursts, we considered a uniform continuum model to assess the homogeneous changes in parameters with the source luminosity.We used absorbed cutoff power law (tbabs×cutoffpl) to describe the continuum emission following the previous studies (Jaisawal et al. 2019;Kong et al. 2020).Following the continuum model, the prominent positive residuals in the energy ranges of 6-7 keV and 0.9-1.1 keV were detected, mainly during the 2017-2018 giant outburst.Moreover, we detected ≈ 7.1 keV iron edge feature only during the high luminosity phases.To account for the positive residuals in the 6-7 keV range, we employed one The 0.7-10 keV energy spectrum of Swift J0243.6+6124obtained from the NICER observation on MJD 58065 (ID 1050390115) near the peak of the X-ray outburst.The second, third, and fourth panels from the top show the evolution of residual after fitting the continuum, and subsequent addition of 1 keV, 6-7 keV Fe-line complex, and edge feature, respectively.
to three Gaussian components depending on the source luminosity, as recommended in Jaisawal et al. (2019).The residuals in the 0.9-1.1 keV band were modeled using a single Gaussian component.Figure 7 shows representative spectral modeling of pulsar energy emission and spectral residuals after the continuum and line emission components in different panels, obtained from a NICER observation near the peak of the giant outburst in 2017 around MJD 58065.
The continuum parameters obtained from the bestfitted models on all NICER observations from multiple outbursts are presented in Figure 8.The variation of spectral parameters such as hydrogen column density (N H ), photon index (P I), and cutoff energy (E cut ), with unabsorbed luminosity, is shown in panels (a), (b), and (c) of Figure 8, respectively.The uncertainties in the parameter values were calculated within the 90% confidence range.Subsequently, the luminosity was estimated from the unabsorbed flux in the 0.7-10.0keV energy range, assuming a distance of 7 kpc (Wilson-Hodge et al. 2018).
Furthermore, we noticed the N H varies in a narrow range between (0.9-1.2) × 10 22 cm −2 across these observations.Therefore, to restrain any spectral degeneracies, we refitted the spectra using an absorbed cutoff power-law model at a fixed average value of N H =1.064× 10 22 cm −2 .Following this approach, a positive residual near the lower energy side arises that could be described with a soft bbodyrad component.The luminosity variations of the latter model parameters obtained from an absorbed cutoff power-law model with a blackbody component are shown in Figure 9.
By examining the overall behavior of these parameters from Figures 8 & 9, two transition points can be identified.The photon index (P I) exhibited a two-component transition with luminosity that may be associated with the changes in the accretion mode (Figure 8).The first transition occurred around a luminosity L 1 of 7.5×10 37 erg s −1 , where P I showed an increasing trend with luminosity.A decrease in the photon index was observed up to a luminosity of 2.1×10 38 erg s −1 (L 2 ) which can be identified as a second transition point.The photon index remains almost stable beyond L 2 .A similar trend can also be seen in Figure 9.The transitional luminosities L 1 and L 2 are marked with dotted lines in these figures.
Moreover, we found the cutoff energy (E cut ) increases with luminosity clearly below the first transition point L 1 (Figures 8 & 9).Due to the limited bandpass of NICER, an upper limit of 30 keV was imposed on the cutoff energy during fitting.This is in line with the maximum value observed in this pulsar based on broadband spectral analysis using HXMT data (Kong et al. 2020).Between L 1 and L 2 , the E cut steeply decreased with luminosity, whereas a gradual evolution is observed above L 2 .Furthermore, we would like to highlight the unique evolution of P I & E cut with luminosity observed beyond the second transition point in our study.These parameters take a distinct trajectory compared to the trends reported in Figure 3 of Kong et al. (2020) in the super-Eddington regime.While Kong et al. (2020) identified a positive correlation between P I & E cut and luminosity, our findings indicate a constant P I and an anti-correlation trend for E cut with luminosity.It is important to note that in the super-Eddington phase, the emission from an outflow or a reflection component can alter the shape of the X-ray continuum.Such an effect may not be tracked alone with NICER due to its limited energy coverage.Therefore, we can expect a slight variation in parameters in the super-Eddington regime.Additionally, to examine this, we conducted spectral fitting by introducing an additional blackbody component during the super-Eddington phase, as suggested by Tao et al. (2019).The temperature of this blackbody component exhibited an evolution between ≈1 and 3 keV, potentially stemming from the contributions of the hotspot and top of the accretion column, with a radius between 10-40 km.In the presence of the second blackbody component, the photon index and cutoff energy show a positive pattern after the second transition point, similar to Kong et al. (2020).However, it is crucial to exercise caution in interpreting these findings due to the constraints imposed by the limited energy band of NICER.
In panels (c) and (d) of Figure 9, the evolution of blackbody temperature and its corresponding emission radii are presented.Below L 2 , the blackbody temperature (kT ) gradually varies in a range of 0.1-0.7 keV.The size of the emission site also changes between 10 to 30 km, for the source distance of 7 kpc, below L 2 .A sudden change in these parameters is observed around L 2 .
In addition to evolving iron emission lines in the 6-7 keV band (Jaisawal et al. 2019), an emission line around 1 keV is detected in the spectra.The 1 keV feature is mainly detected during the giant outburst at luminosities above 2×10 38 erg s −1 i.e. the second transition point.The parameters obtained from the emission line analysis are presented in Figure 10.From this figure, it can be observed that the variation of central energy with luminosity falls within the error bars.However, the line width (σ) and equivalent width (EW ) of the 1 keV line show an increase as the luminosity increases up to ∼8×10 38 erg s −1 , and beyond that decreases with luminosity.This indicates a strong luminosity dependency of the 1 keV emission line.

DISCUSSION
We have studied the first Galactic ultraluminous Xray pulsar Swift J0243.6+6124 using NICER observations covering a giant and multiple normal X-ray outbursts between 2017 and 2023.Such coverage allows us to probe the timing and spectral properties of this unique pulsar depending on the mass accretion rate.Our analysis over multiple outbursts revealed a wide luminosity variation in a range of (0.1-153)×10 37 erg s −1 in 0.7-10 keV energy range, assuming a distance of 7 kpc.The discussion on timing results obtained from normal outbursts is presented first in this section.We thereafter discuss the implications of observed two spectral transitions from the source.L 37 (erg s -1 ) Figure 9.
The panels (a)-(d) show the luminosity dependencies of spectral parameters such as photon index (P I), cutoff energy (Ecut), blackbody temperature (kT ), and blackbody radius at a fixed column density over multiple outbursts of Swift J0243.6+6124.The symbols have the same meaning as Figure 8.
We examined the spin frequency evolution of the pulsar during its recent 2023 outburst (Figure 2).The spin frequency evolution from previous outbursts has been reported by Wilson-Hodge et al. (2018) and Serim et al. (2023).Before the 2023 outburst, the source was in a quiescent state between MJD 58535-60097, where a decrease in spin frequency was observed from 102.10 to 102.03 mHz as per the Fermi/GBM.However, the source gradually spun up to 102.12 mHz after the 2023 outburst (Figure 2).This value is similar to the spin frequency of the neutron star observed after the giant and subsequent normal outbursts between 2017-2019(Wilson-Hodge et al. 2018;Liu et al. 2023;Serim et al. 2023).A significant amount of mass transfer is expected to enhance the spin frequency due to the transfer of angular momentum to the neutron star.
We further have investigated the pulse profile of the pulsar from post-giant outbursts including its recent 2023 outburst (Figures 3).The pulse profiles offer insights into the geometry of the emission regime on the neutron star surface.The pulse profiles of accretionpowered pulsars tend to be simpler and smoother in hard X-rays.However, the soft X-ray pulse profiles appear complex due to the influence of circumstellar scattering and absorption (White et al. 1983).In our analysis, we observed complex pulse profiles in the 0.5-10 keV energy band , featuring various dips or notches at different pulse phases (see Figure 4).Similar features have been observed in the pulse profiles of other Be/X-ray binary pulsars such as V0332+53, 1A 0535+262, EXO 2030+375, GX 304-1, and RX J0209.6-7427(Tsygankov et al. 2006;Naik et al. 2008Naik et al. , 2013;;Naik & Jaisawal 2015;Epili et al. 2017;Jaisawal et al. 2016;Vasilopoulos et al. 2020).These features are usually attributed to the absorption of photons by matter streams locked at specific pulse phases of the neutron star, demonstrating the dynamics of matter distribution in the magnetosphere.
The critical luminosity for Swift J0243.6+6124 can be considered to be ≈10 38 erg s −1 based on the evolution of the source during its giant outburst (Wilson-Hodge et al. 2018).Assuming this limit, the pulsar was accreting in a sub-critical or close to critical luminosity regime during the normal outbursts between MJD 58303-58533 and the recent 2023 outburst where the observed luminosity was in the range of (0.2-9.0)×10 37 erg s −1 .We presented the pulse profile evolution of the pulsar with luminosity during the normal outbursts in Figure 3. Initially, the pulse profiles are single peak dominated in nature when the luminosity was below ∼0.5×10 37 erg s −1 .Between ∼0.5-6×10 37 erg s −1 , various dips/notches arise and the profile evolves to double-peaked.A clear double-peaked profile is visible below the first transition point L 1 .Moreover, the pulse profile evolved from a double-peaked to a smooth singlepeaked structure at luminosities above ∼6×10 37 erg s −1 , which is similar to the pulse profiles shape reported by (Wilson-Hodge et al. 2018) around this luminosity.The observed pulse profile evolution even below the critical regime suggests a change in the emission geometry or beam pattern of the neutron star depending on the mass accretion rate.
In addition to the pulse profile, the pulsed fraction (P F ) provides valuable insights into the pulsating emissions originating from a source.We observed a variation in the P F corresponding to changes in luminosity (Figure 5).The pulsed fraction (P F ) exhibits a moderate variation from ∼ 10% to 15% within the luminosity range of ∼ 2×10 36 to ∼ 9×10 37 erg s −1 .Notably, during the giant outburst, Wilson-Hodge et al. (2018) found an increase in the pulsed fraction with rising luminosity above a critical luminosity of 10 38 erg s −1 from 20% to 55%.Below the critical luminosity, the P F varies between 15-30% during the giant outburst.This is similar to our present findings during the post-giant outbursts.Based on the examination of P F and observed outburst luminosities, the source was accreting below or close to the critical regime at the peak of these multiple normal outbursts.

Detection of QPO and break feature in the PDS
We observed low-frequency QPOs from Swift J0243.6+6124 at particular epochs within a given range of luminosity during the subsequent normal outbursts after the giant outburst.The detected QPO frequencies range from 30 to 47 mHz within a luminosity range of (0.2-2.0) × 10 37 erg s −1 (see Table 1).In panels (b) and (c) of Figure 5, we show the evolution of QPO frequency and width with the luminosity, respectively.The QPO frequency appears to be positively related to luminosity overall, however, the width remains almost flat.We detected no QPO-like signatures in the 2023 X-ray outburst with NICER.No QPO was also found in the NuSTAR observation during the 2023 outburst (Pradhan et al. 2023).During the 2017-18 giant outburst of Swift J0243.6+6124, the QPOs with frequencies between 50-70 mHz within the luminosity range of 0.28-2.05× 10 37 erg s −1 in 0.2-12 keV was observed with NICER by Wilson-Hodge et al. (2018).Moreover, the power spectra obtained with Insight-HXMT in the 20-40 keV range revealed weak QPOs with luminosity-dependent frequencies ranging from ≈50 to 200 mHz during the giant X-ray outburst (Doroshenko et al. 2020).This shows the transient nature and energy dependency of the low-frequency QPOs.The detected QPO frequencies in our study are consistent with the previous QPO studies in other high mass X-ray binary (HMXB) accreting pulsars (Paul & Naik 2011).Also, the phenomenon of QPO detection primarily in the low luminosity phase has been observed in other cases, such as in KS 1947+300 (James et al. 2010).
In the case of HMXBs, the physical origin of these QPOs is commonly explained by two models: the Keplerian Frequency Model (KFM; van der Klis et al. 1987) and the Magnetospheric Beat Frequency Model (MBFM; Alpar & Shaham 1985).The QPO feature arises when the accretion process is governed by the interaction between the co-rotating magnetosphere and the inhomogeneities in the inner accretion disc, which creates variabilities in the mass accretion rate.KFM states this variable mass accretion rate occurs at the Keplerian frequency.For MBFM, this happens at a beat frequency between pulsar spin frequency and mat-ter rotational frequency at the inner disc.In the case of Swift J0243.6+6124, the spin frequency of the pulsar is approximately 102.1 mHz, while the detected QPO frequency ranges from 30 to 47 mHz, nearly one-third of the pulsar frequency.When the spin frequency of the pulsar exceeds the Keplerian frequency at the inner edge of the accretion disc, the co-rotating magnetosphere throws away the accreted matter outwards, which is also called centrifugal inhibition of accretion (Lamb et al. 1973).The KFM is therefore applicable only when the QPO frequency is higher than the neutron star spin frequency, as observed in transient Be/X-ray binary pulsars like EXO 2030+375 (Angelini et al. 1989) and 1A 0535+262 (Finger et al. 1996).Therefore, the MBFM could potentially explain the origin of the quasiperiodic variabilities in Swift J0243.6+6124, which applied in previous studies like 4U 0115+634 (Dugair et al. 2013) and KS 1947+300 (James et al. 2010).
Based on the detected QPO frequency, we can also calculate the magnetic field of the neutron star using the MBFM.The magnetic field of the neutron star for a dipolar structure is given by Where, M 1.4 =1.4M ⊙ (mass of the neutron star), Ṁ−8 =10 −8 M ⊙ yr −1 (mass accretion rate), and R 6 =10 6 cm (radius of neutron star).
According to MBFM, Ω K = Ω QP O +Ω s , where Ω K , Ω QP O , and Ω s are the Keplerian frequency, QPO frequency, and pulsar spin frequency respectively.Within the luminosity range of 0.2-2.0 × 10 37 erg s −1 (assuming a distance of approximately 7 kpc), the mass accretion rate ( Ṁ ) is estimated to be between 0.02-0.13× 10 −8 M ⊙ yr −1 .By applying these values in Equation 2 and considering the detected QPO and pulsar frequencies, the magnetic field (B) of the neutron star is estimated to be approximately 2.5 × 10 12 G.The obtained magnetic field is in good agreement with measurements provided by previous studies using indirect methods (Tsygankov et al. 2018;Doroshenko et al. 2020;Bykov et al. 2022).
For example, the nondetection of the propeller effect by Tsygankov et al. (2018) constrained the upper limit of the magnetic field to 6×10 12 G.Additionally, Doroshenko et al. (2020) provided a range of magnetic field values after analyzing various features such as propeller effect, state transition luminosity from sub-to super-critical state, transitional luminosity from gas pressure dominated (GPD) state to radiation pressure dominated (RPD) state.The calculated magnetic field range was found to be approximately (3-9)× 10 12 G, with the lower limit being more likely.Our magnetic field value is also consistent with the findings of Jaisawal et al. (2019), where the significant broadening of the iron line at higher luminosities indicated the possibility of a lower magnetospheric radius, and hence a relatively low magnetic field ranging from 10 11 to 10 12 G.Bykov et al. (2022) also calculated the magnetic field using the reflection model relxilllp during the super-Eddington phase of the source.They found that during the peak of the outburst, the inner accretion radius was approximately 2-3×10 7 cm, corresponding to a magnetic field of 3×10 12 G.However, some studies have reported a magnetic field of the pulsar as an order of magnitude higher than the estimated value in our work and previous studies.For instance, Kong et al. (2022) detected a cyclotron resonance scattering feature (CRSF) in the 120-146 keV range through phase-resolved spectroscopy during the bright phase of the outburst.The calculated magnetic field corresponding to this observation is approximately 1.6×10 13 G, which is the strongest ever detected for a neutron star in binaries.This component is thought to represent the quadrupole component of the magnetic field.
A noticeable break is also observed in the power density spectra, as shown in Figure 6.The evolution of the obtained break frequency (Br f ) with luminosity is shown in panel (e) of Figure 5. Below the luminosity of ∼7.5×10 37 erg s −1 , the Br f remain around 80 mHz.However, the break frequency reached 140 mHz, exceeding the pulsar spin frequency (101 mHz) beyond a luminosity of 7.5×10 37 erg s −1 .This kind of evolution was also observed by Doroshenko et al. 2020 (Figure 6) in the 20-40 keV Insight-HXMT light curve.The important point to note is that the first spectral transition was also observed around 7.5×10 37 erg s −1 , although the connection is not clear yet.Again, the slope of the power laws below and above the break frequency varies with luminosity between 0.5-1 and 2-3, respectively.Above mentioned slopes Γ 1 & Γ 2 in the PDS are mentioned to be created by variabilities in the accretion disc and magnetosphere of the neutron star, respectively (Hoshino & Takeshima 1993;Revnivtsev et al. 2009).This also suggests the conversion of matter flow from accretion disc flow to magnetospheric flow.The higher values of Γ 2 than Γ 1 indicates suppression of variabilities within the magnetosphere (Revnivtsev et al. 2009).Also, the observed break frequency, being almost 1.4× the spin frequency, suggests that the variability timescale relative to the break frequency is not directly associated with the expected Keplerian timescale at the inner edge of the disc (Revnivtsev et al. 2009).This conclusion can be drawn particularly at low luminosity levels (as in our study), where the spin period is expected to be around the same duration.

Spectral characteristics of Swift J0243.6+6124
To gain a detailed understanding of the X-ray emission mechanisms from the source and to complement our timing studies, we performed spectral analysis using high-cadence NICER data.Our study covers a wide time range from 2017 to 2023 that included the giant outburst, multiple subsequent normal outbursts, and the recent 2023 X-ray outburst.During our spectral analysis, the highest luminosity of the source in the 0.7-10 keV range was estimated to be 1.53×10 39 erg s −1 .This result signifies that during the observation period, the source exceeds the Eddington luminosity limit (L Edd = 1.25×10 38 erg s −1 for a typical 1.4M⊙ neutron star), the maximum possible luminosity for a spherically symmetric emitting source.We are observing luminosity beyond the Eddington limit because of its calculation assumptions: First, the accretion is spherical, and second, the incoming matter and outgoing photons interact through Thomson scattering.But in the case of highly magnetized neutron stars, the interaction cross-section can be much below the Thomson scattering cross-section, and the radiation can escape in the perpendicular direction to the mass accretion direction, which on aggregate can increase the luminosity of the pulsar beyond the Eddington luminosity (Basko & Sunyaev 1976).
Our spectral analysis revealed two distinct transitional luminosities, denoted as L 1 and L 2 , which are ≈7.5 × 10 37 and ≈2.1×10 38 erg s −1 in 0.7-10.0keV energy range, respectively for a distance of 7 kpc.These transitions are observed in continuum parameters such as photon index and cutoff energy, as shown in Figures 8  & 9. Insight-HXMT observed the transition around 1.5 ×10 38 erg s −1 and 4.4 ×10 38 erg s −1 in the 2-250 keV band at a distance of 6.8 kpc (Kong et al. 2020).These two transitions in spectral parameter evolution are observed for the first time in Swift J0243.6+6124.This type of behavior has possibly not been observed before for any other accretion-powered pulsars.In a study conducted by Reig & Nespoli (2013), the spectral analysis of nine Be/X-ray binary pulsars revealed the presence of two branches, indicating a single-luminosity transition.However, the photon index of these sources evolved opposite to Swift J0243.6+6124,i.e., the photon index decreased with an increase in luminosity up to the first transition point (Figure 6 of Reig & Nespoli 2013) except for Swift J1626.6-5156,where the photon index remained constant.This shows that the X-ray emission mechanisms in Swift J0243.6+6124 are distinct from other accretion-powered pulsars.Reig & Nespoli (2013) only found a single transition point in the evolution of spectral parameters with luminosity, representing the presence of two accretion modes on either side of the transition point or critical luminosity.However, in the case of Swift J0243.6+6124,we found two transition points that suggest the presence of three distinct accretion modes, implying more complex behavior in this source.
In accretion-powered X-ray pulsars, the photon index evolves with the luminosity depending on the accretion regimes, where the accretion column acts as the primary source of the X-ray emission.The pulsars usually show a negative correlation between the photon index and luminosity below a single transition point or critical luminosity.A positive correlation is expected between these two parameters in the super-critical luminosity domain.In the sub-critical luminosity regime, the shock region in the column is much closer to the neutron star surface (Becker et al. 2012).The shock region moves up when the mass accretion rate increases, contributing to a taller column in the super-critical luminosity domain.In the latter case, the photons could not acquire enough energy through bulk Comptonization to produce high-energy photons, leading to a softer spectrum with increasing luminosity.On the other hand, in the sub-critical regime (below the transition point), the size of the interaction region decreases with luminosity.This leads to an increase in the optical depth, thereby hardening the spectrum.However, in the case of Swift J0243.6+6124,we made a remarkable observation: there is a positive correlation between the photon index and the X-ray luminosity up to a transition point of approximately 7.5×10 37 erg s −1 (first transition zone).Beyond this point, the correlation becomes negative up to a certain threshold (2.1×10 38 erg s −1 ; second transition zone), and subsequently, we observed no correlation beyond the luminosity of 2.1×10 38 erg s −1 .
First accretion mode (L X ≤ L 1 ): We observe a positive correlation between the photon index and luminosity in the first transition zone below ∼7.5×10 37 erg s −1 .The presence of the observed softer spectrum indicates that lower energy X-ray photons are contributing more to the pulsar energy continuum in comparison to bulk Comptonization of photons with infalling electrons in the accretion column.The softer spectral shape may arise due to blackbody emission from the source.We observed a rise in the soft blackbody (kT ∼ 0.1-0.7 keV) photon emitting region of size (Radius) ∼ 10-30 km in this accretion mode.The observed temperature and size of the emitting region suggest a potential combination of photons from both the neutron star's surface and the hot spot region (Zhao et al. 2019;Elshamouty et al. 2016).This type of soft photon emitting region of temperature between 0.3-0.4keV and size 25-38 km was also identified by Beri et al. 2021 using AstroSat during this luminosity range.Moreover, earlier studies have revealed the presence of a single (Jaisawal et al. 2018;Kong et al. 2020) or multiple (Tao et al. 2019) blackbody components representing the soft X-ray emission from the column and/or optically thick outflow in the spectrum of Swift J0243.6+6124.The presence of these thermal regions may contribute to a softer spectrum.Again as the X-ray luminosity increases, the gas shock region (photon emitting region) is expected to rise within the polar region (panel (b) of Figure 1 in Becker et al. 2012).At this stage, the temperature of the plasma increases in the shocked region.Consequently, the shock height rises with mass accretion rate, resulting in a lower optical depth and less efficient cooling of the accreting plasma.As a consequence, the plasma temperature in the accretion column rises, explaining the correlation between P I and E cut in the first transition zone.
Second accretion mode (L 1 ≤ L X ≤ L 2 ): In this regime, both the photon index and the cutoff energy exhibit a correlated behavior with each other, but both parameters show an anti-correlation with the luminosity.This suggests the presence of a different emission mechanism compared to the previous accretion mode.In the above prescribed luminosity range, the accretion column shows different properties which are also reflected in the timing analysis (Kong et al. 2020).This type of behavior has been previously observed in the sub-critical regime of several sources such as 4U 0115+63, 1A 0535+262, 1A 1118-612, GRO J1008-57 (Reig & Nespoli 2013), EXO 2030+375 (Epili et al. 2017;Jaisawal et al. 2021), 2S 1417-624 (Gupta et al. 2018), and SMC X-2 (Jaisawal et al. 2023).In this luminosity range, it is expected that a radiation shock is present in the accretion column, though its strength is not sufficient enough to bring the matter to a complete rest at the stellar surface (Basko & Sunyaev 1976;Becker et al. 2012).Instead, Coulomb interactions near the base of the accretion column play a significant role in reducing the plasma velocity (Burnard et al. 1991;Nelson et al. 1993).As a result, the height of the emitting region decreases with increasing luminosity following the relation of h e ∝ L −5/7 X (Equation 51 in Becker et al. 2012).The reduction in the size of the sinking region, or the Comptonization region, leads to an increase in optical depth, resulting in the production of harder photons and thus a lower photon index (Becker et al. 2012).Additionally, the cutoff energy decreases as the cooling mechanism through Comptonization dominates over the heating mechanism in this lumi-nosity range, given the increased density resulting from the decreasing height of the emission region.
Third accretion mode (L X ≥ L 2 ): In this luminosity range, the source is in a super-critical state, and we found that the P I remains almost constant with luminosity, while the E cut value decreases.The super-critical state is characterized by the dominance of radiation pressure in controlling the overall flow dynamics of the plasma, ultimately causing the plasma to come to rest on the stellar surface (Davidson 1973;Basko & Sunyaev 1976).The height of the emission region in this case increases with luminosity as h e ∝ L X (Equation 40, Becker et al. 2012).In this regime, the effective velocity of the incoming electrons decreases due to the balanced advection (inwards) and diffusion (outwards) of photons (Becker et al. 2012).As a consequence of the low effective electron velocity, photons do not acquire sufficient energy through bulk Comptonization, leading to spectral softening.Moreover, in this state, we detected a blackbody emission with a temperature of ∼0.1 keV and size of ∼200 km.The blackbody with such characteristics may indicate the presence of optically thick outflow during the super-Eddington phase as suggested by Tao et al. (2019) & Beri et al. (2021).
In our spectral analysis within the energy range of 0.7-10 keV, in addition to the evolving iron emission lines (see e.g., Jaisawal et al. 2019), we also observed the presence of emission lines around 1 keV (Figure 7).These emission lines exhibited a detectable intensity during the super-Eddington phase.There is a low probability that the emission line is of instrument origin as suggested by the instrument team6 .The ∼1 keV line has been also detected in several sources with NICER such as Serpens X-1 (Ludlam et al. 2018), IGR J17062-6143 (Bult et al. 2021), and NGC 300 X-1 (Ng et al. 2022).The emission line displayed a distinctive singlepeak shape, which led us to model it using a Gaussian function.The variations of the key parameters associated with these lines are presented in Figure 10.Several other studies also reported the presence of this line in ULX pulsars ( Kobayashi et al. 2023 & references therein).These line features may be originated from a blend of Fe-L emissions (Gu et al. 2019).The evolution of the 1 keV line parameters with luminosity is presented in Figure 10.The variation of central line energy with luminosity remained within the error bars.However, σ and EW initially increase with luminosity up to the first transitional luminosity L 1 ( ∼ 7.5×10 37 erg s −1 ), and then beyond that decreases with luminos-ity.The line width reaches up to 0.1 keV.As a result, the velocity of the line-emitting material can be approximated to be around 10% of the speed of light based on Doppler broadening.This may suggest the lines might have originated from the accretion disc or from an ultrafast outflow that has been suggested during the super-Eddington phase of the pulsar (van den Eijnden et al. 2019;Jaisawal et al. 2019).

CONCLUSION
In conclusion, our timing and spectral analysis of the first ultraluminous X-ray source in our Galaxy, Swift J0243.6+6124, using NICER data, provided valuable insights into its behavior during different luminosity phases.We observed a luminosity-dependent break in power density spectra, signaling changing accretion dynamics, and identified quasi-periodic oscillations within a specific luminosity range.During the 2023 outburst, the neutron star exhibited a spin-up state and variations in its pulse profile.Spectral analysis revealed two luminosity-dependent transitions at a luminosity of L 1 ≈ 7.5×10 37 erg s −1 and L 2 ≈ 2.1 × 10 38 erg s −1 in continuum parameters, highlighting three distinct accretion modes during the giant outburst.we detected a soft blackbody component (kT ∼ 0.08-0.7 keV), which underwent a discontinuous transition as the source evolved from a sub-Eddington to a super-Eddington state.Notably, during the super-Eddington state, we observed emission lines around 1 keV, indicating X-ray reflection from the accretion disc or outflow material.

3Figure 1 .
Figure1.Swift/BAT monitoring light curve (red solid points) of the pulsar Swift J0243.6+6124 in the 15-50 keV range, between MJD 58010 (2017-09-17) and 60220 (2023-10-03).The light curve is binned by 1 day time frame.The outbursts that are studied using NICER observations are represented with shaded regions.The dark-violet, violet, and green color mark giant, subsequent normal outbursts, and the 2023 normal outburst phase of the source, respectively.

Figure 3 .
Figure 3.The color-coded map of the evolution of the pulse profile of the pulsar with luminosity during the normal outbursts between MJD 58303 and 60190.The pulse profiles are normalized to have values between 0 & 1.Two cycles are shown for clarity.

Figure 4 .
Figure 4. NICER pulse profiles covering a broad range of source luminosity during the normal outbursts between MJD 58303 and 60190.Two cycles are shown for clarity.Here, L37 stands for 10 37 erg s −1 .
figure, the break frequency Br f varied between 50 and 80 mHz below 7.5×10 37 erg s −1 .Beyond this luminosity, the Br f increases to 140 mHz, surpassing the pulsar frequency (101 mHz).Γ 1 and Γ 2 values varied between 0.5-1 and 2-3, respectively, within the studied luminosity range.After fitting the PDS continuum with a broken powerlaw model, we examined the residuals for the presence of potential QPOs.In addition to the neutron star's spin frequency and its harmonic components, we observed a broad hump-like residual below the pulsar spin frequency in a few cases (around 16 IDs).To explore this further, we employed a Lorentzian function (lorentz) as the shape of this hump appeared asymmetric.The Lorentzian function has been widely employed in several QPO studies (see e.g.,Belloni et al. 2002).Figure 6 displays the best-fitted model for the observation ID 1050390170 (MJD 58322) in the PDS.Furthermore, the significance of the QPOs is determined using the method described inBoutelier et al. (2010) based on the Lorentzian fitting.We found that

Figure 5 .Figure 6 .
Figure5.The luminosity evolution of different parameters from PDS analysis during normal outbursts.Panels (a-f) display the values of pulsed fraction (P F ), break frequency (Br f ), QPO frequency (QP O f ), QPO width (QP OW ), slope of power-law before Br f (Γ1), and slope of power-law after Br f (Γ2) as functions of luminosity.The red line represents the best-fitting curve for these values using the spline interpolation method.The shaded area indicates the moving average standard deviation of data points.The purple data points represent values obtained from normal outbursts, except for the 2023 outburst that is shown in black.

Figure 8 .
Figure 8.The panels (a)-(c) show the evolution of the parameters such as column density (NH ), Photon Index (P I), and Cutoff energy (Ecut) with luminosity obtained from multiple outbursts of Swift J0243.6+6124.Panel (d) shows the relationship between the photon Index (P I) and cutoff energy (Ecut).The filled and open markers represent data points obtained from giant and subsequent normal outbursts, respectively.The black points represent the values obtained from the 2023 outburst.

Figure 10 .
Figure 10.The figure illustrates the evolution of Gaussian model parameters of 1.0 keV emission line.The parameters include the central line energy (E1.0), width (σ), and equivalent width (EW ), presented from top to bottom, respectively.

Table 1 .
Observed QPO frequency, width, and its significance with luminosity.The uncertainties are presented at the 1 σ confidence level.