Pulsed Iron Line Emission from the First Galactic Ultraluminous X-Ray Pulsar Swift J0243.6+6124

We report the phase-resolved spectral results of the first Galactic pulsating ultraluminous X-ray source (PULX) Swift J0243.6+6124, modeling its 2017–2018 outburst peak using data collected by the Hard X-ray Modulation Telescope (Insight-HXMT). The broad energy coverage of Insight-HXMT allows us to obtain a more accurate spectral continuum to reduce the coupling of broad iron line profiles with other components. We use three different continuum spectrum models but obtain similar iron line results. For the first time, we detect the pulse characteristics of the broad iron line in a PULX. The variation in the width and intensity of this iron line with σ ∼ 1.2–1.5 keV has a phase offset of about 0.25 from the pulse phase. We suggest that the uneven irradiation of the thick inner disk by the accretion column produces the modulated variation of the broad iron line. In addition, the nonpulsed narrow line is suggested to come from the outer disk region.


INTRODUCTION
Super-Eddington accreting pulsars can power ultraluminous X-ray sources (ULX), which are among the brightest astrophysical sources (above ∼ 10 39 erg s −1 ) in the X-ray sky (for a review see Kaaret et al. 2017).Since Bachetti et al. (2014) first observed X-ray pulses from an ULX (M82 X-2), many pulsating ULX (PULX) sources have been discovered one after another (Fürst et al. 2016;Israel et al. 2017a,b;Carpano et al. 2018;Chandra et al. 2020).A pulsar in an X-ray binary system is able to produce such a high luminosity due to the accretion geometry created by their compact stellar surfaces and complex magnetic fields: an accretion column can be formed near the magnetic poles of the pulsar, and strong radiation escaped from the sides of the column will produce a fan-beam emission.The accretion column illuminates a large surrounding area, allowing us to study the spatial distribution of the surrounding matter.
Many astronomical X-ray sources produce fluorescent iron lines by reprocessing of illumination from the central engine by surrounding materials such as the accretion disk.The emission lines can be skewed and broadened by the combined effects from Doppler effect, gravitational redshift, and Compton Scattering.Therefore, it is possible for us to analyze the ionization degree, kinematics, and geometric distribution of materials around the radiation area through the spectral structure shown by the iron lines, especially in neutron star systems usually with complex accretion structures.In fact, pulsed iron line emissions associated with rotation are detected in some high mass X-ray binaries (HXMBs) such as Cen X-3 (Day et al. 1993), Her X-1 (Choi et al. 1994), LMC X-4 (Shtykovsky et al. 2017), 4U 1538-522 (Hemphill et al. 2014), GX 301-2 (Liu et al. 2018), and V 0332+53 (Bykov et al. 2021), but none of them are PULX.The known PULXs are distant (> 80 kpc, except our target), and the limited number of the pulsed photons is difficult to support the analysis of the emission line.We need higher statistics to resolve the behavior of the iron line in order to probe the geometry at the super-Eddington accretion state of a PULX.
The transient X-ray source Swift J0243.6+6124 was discovered by the Neil Gehrels Swift Observatory in October 2017, when it reached a flux of ∼ 80 mCrab at the rising phase of a giant outburst (Cenko et al. 2017;Kennea et al. 2017).The detection of the pulse period of ∼ 9.86 s and the optical counterpart of a Be star in this source, identified it as a Be X-ray binary pulsar (Kennea et al. 2017;Ge et al. 2017;Jenke & Wilson-Hodge 2017;Kouroubatzakis et al. 2017;Reig et al. 2020).The distance to this Be star measured by Gaia DR2 parallax is 6.8 kpc (Bailer-Jones et al. 2018).Together with the peak flux exceeding 8 Crab, this distantce reveals that the source is the first Galactic PULX.Such a bright source makes a good laboratory for studying the physics of super-Eddington accretion.
In the super-Eddington state of Swift J0243.6+6124, the broad iron line complex has been resolved into components centered at 6.4 keV, 6.67 keV, and 6.98 keV by Jaisawal et al. (2019) using high energy resolution data from NICER and NuSTAR.The putative neutral 6.4 keV line has a fairly large width, with σ ∼ 1.67 keV when approaching the outburst peak.The edge feature at 7.1 keV is obvious at this time.Kong et al. (2020) reported the evolution of the spectral parameters during the outburst with Insight-HXMT observations, where σ reached 1.7 keV at the peak.They also gave phase-resolved spectral parameters at the outburst peak (Kong et al. 2022), but did not analyze the properties of the line emission in details.Bykov et al. (2022) also studied the phase-resolved spectrum using the reflection model with NuSTAR observations, and obtained a pulsed reflection fraction in the super Eddington state.We note that such a broad iron line is difficult to be separated from the continuum component in the spectrum if the energy range is not wide enough, as the continuum model and cutoff power-law plus blackbody for accretion pulsars may be coupled to the broad line.Therefore, we hope that more accurate phase-resolved line parameters can be obtained by conducting detailed spectral analysis of the broad-band Insight-HXMT observations.In section 2, we list the observations and methods of data reduction.In section 3, we present the spectral results for different models.We discuss our results in Section 4 and make conclusions in Section 5.

OBSERVATIONS AND DATA REDUCTION
Insight-HXMT carries three slat-collimated X-ray telescope: the High energy X-ray telescope (HE), the Medium Energy X-ray telescope (ME), and the Low Energy X-ray telescope (LE).Its broad energy band (1−250 keV), large detection area (5100 cm 2 in 20−250 keV for HE) and the small dead time make it a powerful satellite in X-ray spectral studies of bright X-ray sources (Zhang et al. 2019(Zhang et al. , 2020)).
Insight-HXMT made 102 pointing observations of Swift J0243.6+6124 during its 2017−2018 outburst.The Insight-HXMT Data Analysis Software (HXMTDAS) v2.04 with default filters is used to reduce the data.These include screening data with elevation angle (ELV) > 10 • , geometric cutoff rigidity (COR) > 8 GeV, offset for the point position ≤ 0.04 • , and time beyond 300 s to the South Atlantic Anomaly (SAA).We use HXMTDAS tasks hespecgen, mespecgen and lespecgen to generate the spectra.The background spectra are estimated with hebkgmap, mebkgmap and lebkgmap.The backgrounds are estimated via version 2.0.9 of the current standard Insight-HXMT background model (Liao et al. 2020a,b;Guo et al. 2020).And the response matrices are created by herspgen, merspgen and lerspgen tasks.
We select ten exposures (Table 1) near the brightest epoch from 58060 to 58066 (MJD) as used by Kong et al. (2022), because the spectral parameters remain stable for these ten exposures (Kong et al. 2020).The photon arriving times are corrected by solar system barycenter and binary-orbiting modulation (refer to the parameters in Table 2).The parameters of binary orbit are taken from the website of Fermi/GBM Accreting Pulsar Histories1 .The daily spin frequencies taken by the Fermi/GBM allows us to use cubic spline interpolation to fit the frequency with time ν(t).Considering the changing pulse-period, the phase-coherent pulse profiles can be derived by calculated a sequential pulse phase ϕ(t) for corrected event time t as where t 0 is fixed at 58027.499066 (MJD), which is the epoch of the first Fermi/GBM periodicity detection (Sugizaki et al. 2020).We set ϕ 0 = 0.07 to make the minimum region of pulse profile fall within the phase 0.1−0.2.To generate the phase-resolved spectrum, we use hxmtscreen to further filter the photons within a specific phase interval, and then generate the spectrum following procedure mentioned before.Finally, we combine the spectra, the background and the response corresponding to the same phase interval.We also combine the spectra of these ten observations as the phase-averaged spectrum.
The spectrum is fitted by the XSPEC 12.11.1 software package (Arnaud 1996).We use the grppha to group minimum photon counts for HE, ME, LE to be 200, 50, and 30, respectively.The energy bands for LE, ME and HE used for spectral fitting are 1−10 keV, 10−30 keV and 28−100 keV, respectively.There are some calibration residual structures of Si and Ag around 1.6−1.9keV and 19−23 keV, respectively, which have been ignored when making spectral fitting (Li et al. 2020).Data in these two energy regions are neglected in subsequent fittings.We further set the systematic uncertainty to 0.5% in all energy band on the basis of Kong et al. (2022).The error of a single parameter of interest is quoted at the 90% confidence level.

Phase-averaged spectrum
The model consisting of a single blackbody plus a cutoff power-law is a good approximation for the continuum spectrum of an accretion-powered X-ray pulsar (White et al. 1983;Bildsten et al. 1997).We first apply Model-I: constant × tbabs × (bbodyrad1 + bbodyrad2 + gaussian1 + gaussian2 + cutof f pl) × edge to the phase-averaged spectrum in 1−100 keV.When fitting the hump structure in 10−30 keV, the additional bbodyrad2 (kT ∼ 3.5 keV , R bb ∼ 1 km) can reduce χ 2 by 96 and its necessity has also been mentioned in previous works (Kong et al. 2020(Kong et al. , 2022;;Tao et al. 2019).The Tuebingen-Boulder model tbabs (Wilms et al. 2000) is used to account for the absorption due to the interstellar medium (ISM).After fitting by Model-I, a sharp line feature at 6.7 keV with σ < 10 −4 keV and an emission-like feature around 9 keV with σ ∼ 0.18 keV still remain (Figure 1, b), with intensities approximately 0.02 and 0.25 times that of gauss2, respectively.When combining the integrated spectrum of 10 observations, the systematic error dominates the spectral residuals.Compared to other residual distributions, these structures are weak calibrated residuals (Li et al. 2020), which have little effect on the fitting results of the iron lines.We note that there is a soft excess in 1−2 keV.Soft excess features are common in accreting pulsars with X-ray luminosity L X ≳ 10 38 erg s −1 (Hickox et al. 2004).Therefore, we should ignored the 1−2 keV energy range for Model-I.The residuals are shown by the blue dots of 2−100 keV in Figure 1    Careful fitting of the continuum can help to limit the parameters of broad iron line.The broad iron line complex can be fitted by double Gaussian with an absorption edge.It was reported by Jaisawal et al. (2019) near the outburst peak of Swift J0243.6+6124.The broad line with σ ∼ 1.3 keV is strongly coupled with the blackbody component with kT ∼ 1−2 keV in spectrum fitting.Removing the blackbody model causes the width of the broad line to increase to more than 2 keV.By taking into account of the energy range of 1−2 keV, we can thus obtain more accurate blackbody parameter measurements and thus constraining the broad iron line better.For the soft excess in 1−2 keV, we draw on the work of Tao et al. (2019) to add a blackbody component with a kT ∼ 0.1 keV for Model-II: constant × tbabs × (bbodyrad1 + bbodyrad2 + bbodyrad3 + gaussian1 + gaussian2 + cutof f pl) × edge.et al. 2016), is included to address the observed reflection features in Model-IV: constant×tbabs×tbpcf ×(bbodyrad1+ relxilllp + cutof f pl).Relxilllp models the radiation of the accretion disk reflecting the cut-off power law spectra from a lamp post geometry.Similar to Bykov et al. (2022) in analyzing the NuSTAR observations of this source, the spin parameter, a, velocity of lamp-post source, β, and redshift, z, are all set to 0, and the height of the lamp-post source, h, is set to 5 (in gravitational radii).The outer disk radius, R out , is set to 400 (in gravitational radii).We bind the power-law index, Γ, and cutoff energy, E cut , of the incidence spectrum to the primary spectrum.The reflection fraction, f refl , is then fixed at -1.These 4 models are listed in Table 3 and the parameters of phase-averaged spectrum is list in Table 4.We mark the baselines corresponding to the individual ratios on the Y-axis with black dashed lines, as well as the 6.67 keV and 7.1 keV on the X-axis.We note that the structures shown in the residuals around 10 keV, which is the boundary between LE and ME, likely arise from the energy range of the broad iron line not taken into account by this simple continuum modeling process.The pulse profile by HE data is plotted on the right which brighter position represents higher pulse intensity.

Iron line profile
For an overview of the phase-resolved iron line, we adopt the model: constant × tbabs × (bbodyrad1 + cutof f pl) to account for the continuum of the phase-resolved spectrum.The spectral ratios of the data with the fitted continuum model are shown in Figure 2. On the phase-resolved continuum, the hump structure in 20−30 keV is obvious at phase 0.1−0.6.The parts of the absorption feature are present on the 80−100 keV spectrum at phase 0.7−1.0.
The broad iron line complex dominated by lines peaking at 6.67 keV was reported by Jaisawal et al. (2019) at the outburst peak of Swift J0243.6+6124.In their work, the iron line complex consists of a broad 6.4 keV line, a narrow 6.67 keV line, and a narrow 6.98 keV line.In our fitting of the phase-averaged spectrum, the line energy of the broad iron line is approximately at 6.55 keV (Table 4).The iron K-edge feature at 7.1 keV is also obvious in our ratio plot.But we did not detect the 6.98 keV line structures in Insight-HXMT data.As shown in these ten phase-resolved spectral ratio plots, the iron line emission is stronger at pulse-min.The equivalent width reaches 1 keV at pulse-min but only 0.3 keV at pulse-max.In order to quantify these pulsed properties, we still use the same models with free parameters to measure the iron lines.Previously, Kong et al. (2022) performed a detailed phase-resolved analysis of the cyclotron resonance scattering features (CRSF) at the peak of the outburst.We follow their fitting process but focus on the broad iron line in the lower energy range.We thus take the range of 1−200 keV to fit this component with gabs to obtain more accurate and appropriate continuum parameters.The width of gabs is fixed at 20 keV and thawing it has little effect on the continuum parameters.The parameters of gabs are fixed when the fitting only reaches up to 100 keV.

Phase-resolved spectrum
The bbodyrad2 is only needed in the half cycle around the minor peak (phase 0.1−0.6).The broad hump which peaks at ∼ 20−30 keV can be seen in phase 0.1−0.6 from Figure 2. It may be the Compton hump, but we use the blackbody model to fit this continuum structure.For spectra in other phase ranges, adding an extra blackbody component to the model improves the χ 2 very slightly only, but has a strong coupling effect on the iron line and the determination of other continuum component.In the fitting of the phase-resolved spectra, kT 2 reaches 4−6 keV in phase 0.1−0.6.After averaging with spectra beyond phase 0.1−0.6 which does not contain this component, kT 2 of the phase-averaged spectrum is only about 3.5 keV.The values of n H of the tbabs for Model-I, Model-II and Model-III are fixed at 1.0, 1.1 and 0.7 (based on parameters in Table 4), respectively, to better constrain the variation of the intrinsic spectrum.Taking into account the energy resolution, the width of gauss2, σ Fe2 , is difficult to constrain at phase 0.5−0.8 and therefore was fixed at 0.15 keV (the cross in panel e of Figure 4.).
The continuum spectral parameters are shown in Figure 3.For the fitting results of Model-II, the phase-modulated blackbody component with the inferred radius of about 2000−10000 km given by Model-II seems physically implausible.For Model-III, the column density between 1.5 − 2.8 × 10 22 cm −2 corresponds to about 10 11 g variation in accretion material of the size of a neutron star.The energy range of the data is not yet able to determine this soft component, so it is only used to fit the continuum for simplicity.We only list the best-fit results for Model-III in Table 5, as the physical interpretation of this soft component is beyond the scope of this work.Some of the phase-resolved spectrum fitting parameters of interest by Model-III are shown in Figure 5.
The best-fit iron line parameters are shown in Figure 4. Model-I, Model-II, and Model-III mainly use different components to handle the fitting in 1−2 keV, but all use the spectral model consisting of two Gaussians with an absorption edge to fit for the line complex in the Fe K band.Similar results about the iron line complex are obtained, which means that our results on the iron line emission do not depend on the continuum emission models and are thus reliable.The broad line Fe 1 varies periodically while the narrow line Fe 2 is almost constant.Therefore, we further fixed E Fe2 and σ Fe2 to 6.63 keV and 0.15 keV respectively, in order to obtain a more significant variation feature of Fe 1 .These parameters are shown in Figure 7.
From the fitting parameters of the iron line, it can be seen that E Fe1 varies around 6.4−6.6 keV, but with a large error (Figure 7, a).This is mainly due to the asymmetry of the iron line profile (Figure 2).The broad red wing of the iron line complex will cause the center energy of the broad line to be lower (Jaisawal et al. 2019).Therefore, the line energy and width of Fe 1 appear to be anti-correlated (Figure 4 and 7).Fixing E Fe1 will result in the same variation trends of iron lines, where σ Fe1 is in the range of 1.22 to 1.42 keV, and N Fe1 is in the range of 0.67 to 1.07 (Model-III).The σ Fe1 displays almost a half-period sinusoidal variation and is inversely correlated with the pulse intensity (Figure 7, b).The N Fe1 (count rate of the Fe 1 line) curve leads the pulse profile by nearly 0.25 phase (Figure 7, c).The absorption edge depth, τ (Figure 7, e), presents the same variation trend with σ Fe1 even with a fixed n H (Model-I).We notice the unavoidable coupling between the broad iron line and the iron K-edge.But the maximum improvement in χ 2 , ∆χ 2 , caused by adding an edge component can reach 107 (with 1336 d.o.f. for Model-III) in phase 0.1−0.2, and lower than 20 at the peaks.Changes in the parameter of the edge optical depth, τ , is responsible for this variation.It approaches to 0 in some phase intervals and cause a small value of ∆χ 2 .In order to further study the detailed features of the iron line, we apply reflection model to provide a more physical description of the line profile and continuum.

Reflection modeling of the phase-resolved spectra
The Iron Abundance (in Solar Units) in Model-IV is fixed at 3 to fit the phase-resolved spectra.Phase-resolved parameters of the reflection component in Model-IV are shown in Figure 6, and the other continuum parameters are shown in Figure 3.The complex line profile modeled by the reflection of the accretion disk is a good approximation of the spectrum.There may be additional reprocessing of the primary radiation for relxilllp, which is especially evident in phase range 0.4−0.5.The bbodyrad2 component used here has a maximum flux (Figure 3, panel i and j), and the ionization parameter log ξ is also peaks in this phase range.The parameter log ξ reached 3.6 (Figure 6, panel c) at such a high state.The inner disk radius R in ∼ 70 r g (Figure 6, panel b) is consistent with that measured based on NuSTAR observations (Bykov et al. 2022).The gravitational radius r g = 2GM/c 2 , and r g = 4.13 km for a neutron star with M = 1.4 M ⊙ .The application of reflection model of lamppost geometry provides a plausible explanation for the origin of iron line on accretion disks.In the following we will discuss the possible origin of the iron lines by combining the quantitative estimation from the gauss model and the spectral characteristics of the reflection model.

DISCUSSION
Swift J0243.6+6124reached a luminosity of ∼ 10 L Edd (L Edd = 1.8×10 38 erg s −1 for a neutron star with M = 1.4 M ⊙ ) at its 2017−2018 outburst peak.In this work, we performed a phase-resolved spectral analysis during this brightest period from broad-band Insight-HXMT observations, with the purpose of highlighting the physical properties of this PULX during the extremely high accreting states.The results show a pulsed broad iron line feature whose properties are not sensitive to the continuum spectral models.In the following we use the iron line emission features to constrain the accretion geometry of this PULX, under the framework of the radiation beam (Inoue 2020) from accretion column and compact relativistic thick disk at the ultra-luminous state (Kong et al. 2020;Doroshenko et al. 2020;Bykov et al. 2022).

The broad iron line profile
The detection of iron lines by Insight-HXMT shows that it consists of a broad line centered at around 6.5 keV and a narrow line centered at around 6.63 keV (consistent with Jaisawal et al., but the 6.98 keV line was not detected by Insight-HXMT).Based on the measurements of relativistic reflection model with Insight-HXMT data, the iron lines may originate from the accretion disk with R in ∼ 70 r g and log ξ ∼ 3.63.From this perspective, the broad line could be generated in the dense inner disk region, and the narrow line should be farther away.
The 7.1 keV absorption edge might indicate the presence of neutral iron in optically thin material around the pulsar.It is possible that the requirement of the edge model component in the spectral modeling procedure arises from using the over-simplified gauss model to describe the broad iron line profiles produced at the vicinity of the compact object, which usually have a slow rise and fast decay shape.A detailed discussion regarding the edge model parameters would therefore be an over-interpretation.
The inner radius of 70 r g is consistent with the results obtained by NuSTAR (Bykov et al. 2022).Assuming that the inner disk radius R in is approximately the radius of the magnetosphere, the estimated dipole magnetic field is B = 1 × 10 12 G with a coupling parameter Λ = 0.5 (Mushtukov et al. 2017).The estimated magnetic field is an order of magnitude lower than B ∼ 1.6 × 10 13 G reported by Kong et al. (2022) for the same observations, but with the detected cyclotron resonant scattering feature at around 146 keV.These two magnetic strengths may correspond to dipole and multipole fields, respectively, as also discussed by Kong et al. (2022).However, the accretion disk is most likely truncated by the magnetosphere dominated by a dipole field away from the surface of the neutron star.On the other hand, the detected cyclotron resonant scattering feature should be produced very close to the surface of the neutron star, where the multipole field dominates.
The result of Model-IV is based on the reflection geometry of the lamp post.We chose to use this model because it includes proper treatments of the general relativity effects and the atomic physics relevant to the generation of broad iron lines in accretion disks.However, this model assumes that the accretion disk is symmetrically illuminated by a lamp post source which is probably not the suitable description for the geometry of the accretion columns in neutron stars.For fan beam radiation that is unevenly distributed with phase, the variations of the direct radiation and the reflected radiation are quite out of sync, or even inversely correlated in extreme cases.Deviations of the illuminating source from model assumptions used as inputs in the reflection model may also lead to biased continuum parameter measurements.The occlusion of matter in the co-rotating magnetosphere (accretion curtain) is also not considered in the reflection model.Therefore, we further discuss the origin of the iron line from the periodic modulation of the Gaussian line, so as to avoid the limitation and bias of the Model-IV in describing the different phases of pulsars.

Modulated iron line feature
To discuss the possible radiation distribution of the accretion disk, we use the beam pattern of the polar cone region described by Inoue (2020).When the angle θ R = 80 • (angle of the magnetic axis to the rotational axis) and the inclination of the rotation axis i = 50 • , the simulated pulse profile (Figure 12, b in the work of Inoue 2020) is similar to the structure of main and minor peaks of Swift J0243.6+6124.The radiation of fan beam directed towards the neutron star due to the electron scattering effect and gravitational bending (Inoue 2020;Mushtukov & Tsygankov 2023).Therefore we only draw the main radiation distribution of the fan beam (the yellow light on the inner disk in Figure 7, f and g).Assuming that the accretion column formed at the north and south poles is symmetric with respect to the neutron star, the accretion column illuminates the disk mainly around the same azimuth direction.
Whether the iron line is modulated by the rotation of the pulsar is an important basis to determine its origin.The ratio of the iron line to the continuum (shown in Figure 2) is inversely correlated to the pulse intensity of the continuum emission.This means that the modulation amplitude of the iron line is weaker than the continuum.The iron lines should originate in a wider region than the radiation beam, such as accretion disk.The model proposed by Bykov et al. (2022) provides an explanation of the reflection feature in the accretion disk of this source.
Now we discuss our results on the fluorescent iron line within the framework of this model from observer perspective, as shown in Figure 7. Radially, the broad line is generated mainly by the gravitational red-shift of the ionized line (rather than the neutral line) from the inner disk, while the narrow line is from a more outlying region of the disk.Toroidally, the fluorescent emission of the outer disk is almost constant in different phases, resulting in a non-pulsed narrow line (Figure 4, c).The thick inner discs can be divide into 4 regions in 90 • steps.We named them regions A, B, C, and D and placed region A on the side close to the line of sight of the observer as shown in Figure 7.The regions B and D generate the mainly broadened iron line due to the stronger Doppler effect.Therefore, for the observer, the iron line properties corresponding to the geometry in ϕ 1 (phase 0.6−0.7) is broader than in ϕ 2 (phase 0.85−0.95), as shown in Figure 7.If the accretion column is in ϕ 1 + 0.5 or ϕ 2 + 0.5, it will irradiate regions symmetrical with those when the accretion column is in ϕ 1 or ϕ 2 , resulting in the same width of the iron line.This gives rise to bimodal variation of σ Fe1 .When we consider the thickness of the disk, the fluorescence surface is no longer exposed to the observer uniformly.The inner edge of region A in Figure 7 is obscured by the accretion disk itself, making the iron line in ϕ 2 weaker than in ϕ 1 .However, parameter N Fe1 in ϕ 1 has a large deviation from symmetrical configuration in ϕ 1 + 0.5.This indicates the presence of additional complex accretion structures.
The pulsed blackbody component with a temperature kT of 1 keV may correspond to an optically thick structure in the magnetosphere such as the accretion curtain or the column top (Mushtukov et al. 2017;Tao et al. 2019).This structure could weaken the observed iron line by blocking the direct radiation from the column to the disk and/or fluorescence emission from the inner disk.We notice that in Figure 3 (h) the area of bb 1 in ϕ 1 is quite different from ϕ 1 + 0.5, and thus suspect that the observed blackbody component is mainly around region A (in Figure 7), while the blocking occurs mainly near region C.This might be the reason that the iron line in ϕ 1 is stronger than that in phase ϕ 1 + 0.5.The occlusion of the accretion column top or accretion curtain could also modulate the iron line intensity, but the distribution pattern of material on the magnetosphere is unclear, and so we did not draw this structure in Figure 7.
We suggest that the variation of the Gaussian component is mainly caused by the inhomogeneous illumination and reflection from the thick inner disk as the pulsar rotates.The ionization state should be the same at a certain disk radius.Thus we believe the variation of E Fe1 is most likely a result of the skewed line profile, with the strength of the red wing relative to the blue wing varying with phases (Jaisawal et al. 2019).This may also lead to the pulsed iron K-edge, which is significant as shown by the reduction of χ 2 (138) when adding an edge model in the spectral fitting.Although the similar pulse variations were detected in V 0332+53 (Bykov et al. 2021), also a Be/X-ray transient pulsar, the accretion curtain origin in their work seems out of place here.The reprocessing of optically thick envelope is difficult to form the edge structure, and in such a high state, it is difficult to observe the neutral iron atoms in this (f) The radiation configuration in ϕ1.The radiation from accretion column is mainly concentrated in regions B and D. The apparent velocity is greater in these two regions that produce broader iron lines.(g) The radiation configuration in ϕ2.The radiation from accretion column is mainly concentrated in regions A and C. The apparent velocity is lower in these two regions that make a narrower line width.When considering the thickness of the disk, the fluorescent lines at the inner edge of region A will be obscured, resulting in a weaker iron line.
interior region.Using two Gaussian components and absorption edge can only provide a very rough and most basic estimate of the changes in the line profile.In this paper, we present clearly the phased-resolved variation behavior of the broad iron line in Swift J0243.6+6124 with high S/N spectra, and propose a basic picture to explain the origin of the line variations.However, we will not delve into the detailed physical model in this paper and leave it to future theoretical works.

SUMMARY
We have performed phase-resolved spectral analysis of the ultra-luminous X-ray source Swift J0243.6+6124 in its brightest state with the data observed by Insight-HXMT during the 2017−2018 outburst.In this state, it is characterized by a broad iron line peaked at ∼ 6.67 keV.This profile can be approximated by a narrow Gaussian line centered at around 6.63 keV and a broad Gaussian line centered at around 6.5 keV.Referring to the reflection model, we speculate that the broad line originates from inner disk region (∼ 70 r g ), while the narrow line originates from the outer region.
For the first time, we detect pulsed broad iron line signature in a PULX.The variation in width and intensity of this iron line has a phase offset of about 0.25 from the pulse profile.The variation of narrow line is not significant which supports that it is from the outer disk with almost no modulation.This result is generally consistent with the radiation pattern (Inoue 2020) of the dipole field and the thick disk geometry (Kong et al. 2020;Doroshenko et al. 2020;Bykov et al. 2022).The modulation in the broadening of the iron line is caused by the Doppler effect on the inner disk that is irradiated unevenly in different phases.Modulation in strength of broad iron line may be caused by the occlusion of the thick inner disk itself.We also propose a possibility for modulation caused by occlusion of magnetospheric material.For a more specific speculation on the modulation of the broad iron line and the results of the pulsed iron K-edge, a more detailed model description is still needed.11(1338) (b).The χ 2 fitted by Model-I in the energy range of 1−100 keV is 1850 (1334 d.o.f) while in 2−100 keV is 1371 (1252 d.o.f).

Figure 1 .
Figure 1.The phase-averaged spectrum (a) and reduce residuals fitted with Model-I (b), Model-II (c) Model-III (d) and Model-IV (e).The energy range of Model-I fitting is 2-100 keV, but the residuals is shown in 1-100 keV.The dotted lines represent the iron lines.The solid red line shows the reflection component.The two dashed lines represent the lower and higher blackbody components while the dotted dash line draws the cut-off power-law component.

Figure 2 .
Figure 2. Ratio of observation data to the continuum model (shifts up by 0.2 in turn) in different phase intervals.The continuum is obtained by fitting 2−4, and 9−100 keV spectrum with model: constant × tbabs × (bbodyrad + cutof f pl).We mark the baselines corresponding to the individual ratios on the Y-axis with black dashed lines, as well as the 6.67 keV and 7.1 keV on the X-axis.We note that the structures shown in the residuals around 10 keV, which is the boundary between LE and ME, likely arise from the energy range of the broad iron line not taken into account by this simple continuum modeling process.The pulse profile by HE data is plotted on the right which brighter position represents higher pulse intensity.

Figure 3 .
Figure 3.Comparison of continuum spectral parameters of four models.The green, blue, red, and pink data points represent models I, II, III, and IV, respectively.The pulse profile (black line) is plotted by HE data.(a)-(b) The variation of tbpcf parameters in models III and IV.The column density of equivalent hydrogen, nH, in this model is positively correlated with the pulse profile.The covering factor, fcover, is also positively correlated with the pulse profile.For comparison, the fixed nH of tbabs in models I and II is drawn on panel (a).(c)-(e) The variation of cutof f pl parameters.The photon index,Γ is smaller at the main peak.Both the cutoff energy, Ecut and normalization, Ncut, are positively correlated with the pulse intensity.(f) The goodness-of-fit for these four models.The error is set to 0. (g)-(l) The Blackbody parameters variation.bb1: Both temperature (∼ 1 keV) and radius (∼ 20 km) of the blackbody are positively correlated with the pulse intensity.bb2: The temperature ∼ 5 keV, the radius is maximum in phase 0.4−0.5 (∼ 1.5 km), and this component is only added in model I, II, and III in phase 0.1−0.6.bb3: Temperature (∼ 0.1 keV) is inversely correlated to the pulse, but the radius (∼ 5000 km) is positively correlated.

Figure 4 .
Figure 4. Iron line parameters versus pulsar rotation phase.The parameters obtained with Model-I , Model-II, and Model-III are in green, blue, and red respectively.For reference the black curve represents the pulse profile obtained with the HE data.(a)-(c) Parameters of broad iron line versus rotation phase.(d)-(f) Parameters of narrow iron line.σFe 2 is fixed at 0.15 keV during phase 0.5−0.8 which are marked by cross in panel e. (g)-(h) Parameters of absorption edge.

Figure 5 .
Figure 5.The phase-resolved spectrum and reduce residuals fitted by Model-III (left panel) and Model-IV (right panel).The color of each spectra corresponds to the color of the phase interval in the pulse profile inset.Blue, orange, green, and red represent phase intervals of 0.4−0.5, 0.5−0.6,0.6−0.7,0.7−0.8,respectively.The dotted line draws the iron line (left panel) and reflection component (right panel).The dashed line draws the blackbody component.

Figure 6 .
Figure 6.Parameters of reflection model versus the rotation phase.(a) Inclination angle to the normal of the disk.(b) Inner radius of the accretion disk in gravitational radii.(c) Ionization of the disk in logarithm.The value of 0 for neutral, while 4.7 for heavily ionized.(d) Normalization of the reflection component.The pulse profile (black line) is plotted by HE data.In panels (c) and (d), the parameters in phase 0.4−0.5 have a large shift from the mean value.
Geometric sketch of the broad iron line emission region on the accretion disk of two specific phase intervals from the observer's perspective.The accretion disk is divided equally into region A, B, C, and D toroidally, and let the region A be on the side close to the line of sight of the observer.(a)-(e) The parameters of iron line.The results is fitting with fixed EFe 2 and σFe 2 compared to the parameters in Figure 4. We mark two phase intervals ϕ1 (phase 0.6−0.7)and ϕ2 (phase 0.85−0.95) by the dashed boxes whose radiation configurations are illustrated in panels (f) and (g).

Table 1 .
Information of 10 Insight-HXMT observations used in this work

Table 3 .
Model definitions in this workModel

Table 4 .
Phase-averaged spectrum parameters (García et al. 2014; are shown by the orange dots in Figure1 (c).Adding this model resulted in a reduction of 379 in χ 2 with 1332 d.o.f..If 1−2 keV is ignored, adding bbodyrad3 only results in ∆χ 2 = 13 with 1252 d.o.f.Considering that the accretion materials can block some emission regions, incomplete coverage of the absorption component causes soft photons to escape, which can also produce the soft excess.So we add the model, tbpcf , for Model-III: constant × tbabs × tbpcf × (bbodyrad1 + bbodyrad2 + gaussian1 + gaussian2 + cutof f pl) × edge.The fitted residuals are shown in the green dots in Figure1 (d).Adding this model resulted in a reduction of 339 in χ 2 with 1334 d.o.f..If 1−2 keV is ignored, adding tbpcf only results in ∆χ 2 = 12 with 1252 d.o.f..For comparison with known physical models, the relativistic reflection model, relxilllp(García et al. 2014; Dauser