Phase-dependent Spectral Shape Changes in the Ultraluminous X-Ray Pulsar NGC 5907 ULX1

The discovery of coherent pulsations from several ultraluminous X-ray pulsars (ULXPs) has provided direct evidence of a supercritical accretion flow. However, the geometrical structure of such an accretion flow onto the central neutron star remains poorly understood. NGC 5907 ULX1 is one of the most luminous ULXPs, with a luminosity exceeding 1041 erg s−1. Here we present a broadband X-ray study of this ULXP using the data from simultaneous observations with XMM-Newton and NuSTAR conducted in 2014 July. The phase-resolved spectra are well reproduced by a model consisting of a multicolor disk blackbody emission with a temperature gradient of p = 0.5 (T ∝ r −p ) and a power law with an exponential cutoff. The disk component is phase-invariant and has an innermost temperature of ∼ 0.3 keV. Its normalization suggests a relatively low inclination angle of the disk, in contrast to the previous claims in other literature. The power-law component, attributed to the emission from the accretion flow inside the magnetosphere of the neutron star, indicates phase-dependent spectral shape changes; the spectrum is slightly harder in the pre-peak phase than in the post-peak phase. This implies that the magnetosphere has an asymmetric geometry around the magnetic axis and that hotter regions close to the magnetic pole become visible before the pulse peak.


INTRODUCTION
Ultraluminous X-ray sources (ULXs: Makishima et al. 2000) are point-like sources with their apparent luminosity of L X > 10 39 erg s −1 , often found in the off-nuclear regions of nearby galaxies (e.g., King et al. 2023;Pinto & Walton 2023 for recent reviews).Since the luminosity exceeds the Eddington limit for a typical stellar-mass black hole, ULXs have once been considered to be candidates for intermediate mass black holes (e.g., Colbert & Mushotzky 1999;Zampieri & Roberts 2009).To date, more than 1800 ULXs have been discovered (Walton et al. 2021).
The situation has dramatically changed by the discovery of coherent pulsations from the ULX M82 X-2, indicating that the system consists of a neutron star accret-ing at the super-Eddington rate (Bachetti et al. 2014).Following this discovery, similar pulsations have been detected from other ULXs, e.g., NGC 7793 P13 (Fürst et al. 2016;Israel et al. 2017a), NGC 5907 ULX1 (Israel et al. 2017b), NGC 300 ULX-1 (Carpano et al. 2018), NGC 1313 X-2 (Sathyaprakash et al. 2019).There are approximately 10 ultraluminous X-ray pulsars (ULXPs) that have been identified to date.Since the mass of a neutron star is at most ∼ 3M ⊙ , the extremely high luminosity observed in the ULXPs provides direct evidence for the accretion rate exceeding the Eddington limit, offering an ideal site to study the poorly-understood nature of super-critical accretion flow.
NGC 5907 ULX1 is one of the most luminous ULXPs, located on the disk of the edge-on spiral galaxy NGC 5907, the distance to which is 17.1 Mpc (Tully et al. 2016).The source was first reported in the ULX catalog presented by Walton et al. (2011).Observing this ULX with XMM-Newton and NuSTAR, Israel et al. (2017b) discovered a coherent pulsations with the peak luminosity of ∼ 10 41 erg s −1 .The spin period and its first derivative obtained from the July 2014 observations were P ∼ 1.137 s and Ṗ ∼ −5 × 10 −9 s s −1 , respectively.They also revealed an energy-dependent phase shift in its pulse profiles, which indicates the presence of spectral variations as a function of the pulse phase.However, no quantitative report on the phase-dependent spectral parameters was given in their work.
The spectra of NGC 5907 ULX1 have been revisited in some other literature.Walton et al. (2018a) performed phase-resolved spectroscopy, but the spectral shape of the pulsation component was assumed to be invariable with respect to the pulse phase.In other words, the energy-dependent phase shift reported by Israel et al. (2017b) was not taken into account in their work.Fürst et al. (2017) analyzed the phase-averaged spectra of this ULXP to investigate the properties of the accretion disk.From the relationship between the modeled bolometric luminosity and the observed flux of the disk component, they estimated the inclination angle of the disk to be ∼ 89 • , implying a nearly edge-on view of the system.This result is puzzling, if true, because in a super-critical accretion system, the accretion disk is thought to become geometrically thick to block the line of sight to the pulsating accretion flow inside the magnetosphere when viewed from the side.It should be noted, however, that in spectral analysis of typical ULXP, the spectral parameters of the disk and pulsation components are easily coupled with each other.Therefore, careful spectral modeling is essential to constrain the physical properties of both components.
Motivated by the previous work, we revisit the X-ray spectra of NGC 5907 ULX1 obtained from the July 2014 observations simultaneously conducted by XMM-Newton and NuSTAR.We perform phaseresolved broadband spectroscopy by taking into account the phase-dependent spectral variations (Israel et al. 2017b) to gain insights into the nature of the accretion flow.
The rest of this paper is structured as follows.Details of observations and data reduction are described in Section 2. We present analysis and results in Section 3 and discuss their implications in Section 4. Finally, we conclude this study in Section 5.The errors quoted in the text and table and error bars given in figures represent a 1σ confidence level unless otherwise stated.

OBSERVATIONS AND DATA REDUCTION
In this work, we analyze the archival data of NGC 5907 ULX1 from the European Photon Imaging Camera (EPIC) aboard XMM-Newton (Jansen et al. 2001) and the Focal Plane Module (FPM) aboard NuSTAR (Harri-son et al. 2013).Both instruments have a large effective area and adequate temporal resolution, allowing us to perform high-quality phase-resolved spectroscopy in the wide X-ray band.

XMM-Newton
NGC 5907 ULX1 was observed by XMM-Newton on 2014 July 9 with an effective exposure of ∼ 38 ks for the EPIC-pn (Strüder et al. 2001) and EPIC-MOS (Turner et al. 2001) detectors (ObsID: 0729561301).Both instruments were operating in the Full Frame mode.Since the timing resolution of the EPIC-pn is much better than that of the EPIC-MOS (73.4 ms and ∼ 1 s, respectively), we use only the former in the following analysis.We reprocessed the data using the epchain task in the XMM-Newton Science Analysis System (SAS) v20.0.0.The cleaned event data were corrected to the barycenter of the solar system using the SAS tool barycen by assuming the solar system ephemeris DE200 (Standish 1990).We use only single-and double-pixel events.No significant flare was detected during the observation.Following Israel et al. (2017b), we adopted a source coordinate of RA = 15 h 15 m 58 s .62 and DEC = +56 • 18 ′ 10 ′′ .3, and extracted the source data from a 30 ′′ -radius circular region around the source position.The background data were extracted from an annulus region with inner-and outer-radii of 30 ′′ and 90 ′′ , respectively.

NuSTAR
NuSTAR observed the source with an exposure of ∼ 57 ks using the two X-ray detectors FPMA and FPMB (ObsID: 80001042002), simultaneously with the XMM-Newton observation.The raw data were reprocessed using the NuSTAR Data Analysis Software (NuS-TARDAS, v2.1.2included in HEASOFT V.6.31.1) and the calibration database (CALDB) version 20230420.We used the nuproducts task to generate barycentercorrected light curves and spectra.We extracted the source data from a circular region with a radius of 49 ′′ and the background from an annulus with an outer radius of 98 ′′ surrounding the source region.

Timing Analysis
We perform a pulsation search using the EPIC-pn, FPMA, and FPMB light curves.Only 3-10 keV events are used for this analysis, because all the three detectors have sensitivity to this energy band.We apply an accelerated Fourier method, which searches a grid in the frequency f and frequency-derivative ḟ space.To do this, we introduce the HENaccelsearch task in the HENDRICS software package (Bachetti 2015).Subsequently, the best combination of (f, ḟ ) obtained from this search (for each detector) is analyzed in more detail using the Z 2 n statistics (Buccheri et al. 1983) available in the HENzsearch task.The resulting period and its derivative are as follows: P = 1.137 43(2) s and Ṗ = −5.38 +1.02 −0.92 × 10 −9 s s −1 at 56847.9 (MJD) for the EPICpn, P = 1.137 66(1) s and Ṗ = −5.18+0.18  −0.17 × 10 −9 s s −1 at 56847.4 (MJD) for the FPMA, and P = 1.137 66(1) s and Ṗ = −5.20 +0.18  −0.16 × 10 −9 s s −1 at 56847.4 (MJD) for the FPMB.The obtained values of P and Ṗ are consistent among the three detectors when the reference epoch is aligned, and also agree with the previously reported values of Israel et al. (2017b).Since the most stringent constraints on P and Ṗ are obtained for the FPMA, we use their mean values to fold the light curves of all the three instruments.
Figure 1 shows the resulting folded light curves in different energies, obtained from the (a) EPIC-pn, (b) FPMA, and (c) FPMB.In the profiles for the whole energy range (i.e., 0.2-12 keV for the EPIC-pn and 3-30 keV for the FPMA/FPMB), the peak and minimum of the pulse are found at phase 0.4-0.5 and phase 0.9-1.1,respectively.The energy-dependent phase shift is clearly seen in the EPIC-pn profiles, confirming the previous results of Israel et al. (2017b).We also confirm the larger pulse fraction achieved at higher energies, suggesting that the hard X-ray emission predominantly originates from the accretion flow formed along the magnetic fields of the neutron star.

Spectral Analysis
To perform phase-resolved spectroscopy, we divide the pulse period into 10 phase bins with equal width, and extract spectra from each bin.We then merge the spectra from phases 0.1-0.2 and 0.2-0.3 to represent the prepeak part and phases 0.6-0.7 and 0.7-0.8 to represent the post-peak part of the pulse cycle.The luminosities at these two phases (0.1-0.3 and 0.6-0.8)are nearly identical to each other.We use the EPIC-pn data in the 0.3-10 keV band and the FPMA/FPMB data in the 3-30 keV band for the spectral analysis.The unfolded spectra of the pre-peak and post-peak phases are shown in the top panel of Figure 2 (red and blue, respectively), confirming the harder spectrum in the former, as implied in the pulse profiles (Figure 1).In Figure 6 in Appendix A, we also show the spectra extracted from the eight different phases (i.e., the pre-peak, post-peak, and other 6 phases with the 0.1 phase width).
For more quantitative studies, we fit the spectra with a model consisting of a diskpbb component and an power law with an exponential cutoff (cutoffpl) using the XSPEC analysis package version v12.13.0c (Arnaud 1996).The diskpbb model represents multicolor blackbody emission from an optically-thick accretion disk with a radial temperature gradient (Mineshige et al. 1994).The free parameters of this component are the innermost temperature (T in ), the index of the radial dependence of the disk temperature (p, where T (r) ∝ r −p ), and the normalization (N ).Theoretically, the index p E fold (keV) a 4.9 +0.9 −0.7 4.6 +0.9 −0.7 4.9 +1.0 Notea Uncertainties are obtained using XSPEC steppar command in the 2D-spaces between NH and N , and between Γ and E fold .
b Flux in the 0.3-40 keV.
is expected to range from 0.5 to 0.75, depending on the mass accretion rate and the resulting physical properties of the disk.When the accretion rate is below the Eddington limit, a geometrically-thin, standard disk is formed due to efficient radiative cooling, and p = 0.75 is expected accordingly (Shakura & Sunyaev 1973).When the accretion rate is high, on the other hand, the disk increases its geometrical thickness due to strong radiation pressure and advection, leading to a flatter temperature gradient (p = 0.5) achieved within the radius called "spherization radius" (e.g., Poutanen et al. 2007).This type of disk is commonly referred to as "slim disk" (Jaroszynski et al. 1980;Abramowicz et al. 1988;Watarai et al. 2000Watarai et al. , 2001)).The other model component, cutoffpl, is responsible for the pulsed emission associated with the accretion flow rotating with the magnetosphere of the neutron star.This phenomenological model is often introduced to fit typical ULXP spectra (e.g., Brightman et al. 2016;Walton et al. 2018b), although the radiation mechanism related to this component is still under debate (see §4.3).The free parameters are the photon index (Γ), the e-folding energy of exponential cutoff (E fold ), and flux in the 0.3-40 keV band (F ).
In addition, we introduce the Tuebingen-Boulder interstellar medium absorption model, tbabs, to account for the foreground absorption, assuming the solar abundances of Wilms et al. (2000).The hydrogen column density (N H ) for the Galactic absorption is fixed to 1.38 × 10 20 cm −2 (Kalberla et al. 2005), whereas the value for the absorption in the NGC 5907 galaxy is left as a free parameter.We assume that the absorption intrinsic to the source binary system is negligibly small compared to that in its host galaxy, since this galaxy is edge-on.Lastly, we introduce a cross-normalization parameter, const, to account for the systematic uncertainty in the effective area calibration among the instruments.The const value for the EPIC-pn (C pn ) is fixed to unity, whereas those for the FPMA (C FPMA ) and FPMB (C FPMB ) are allowed to vary freely.
With the model described above, we fit the eight phase-resolved spectra given in Appendix A simultaneously, assuming Poisson data with Gaussian background (PG-statistic : Cash 1979).Since the emission from the accretion disk is considered to be phase-invariant, the parameters of the diskpbb component is tied among the eight spectra.The parameters of the cutoffpl components, responsible for the pulsed emission, are fitted independently.We then obtain a reasonably good fit with the P G-stat value of 7426.11(for 9415 dof), but the parameter p is not constrained within the allowable range of 0.5-0.75.Therefore, we repeat the fitting by fixing the p value to 0.5 or 0.75.The two extreme cases result in similarly good fits with P G-stat values of 7426.10 and 7426.43 for p = 0.5 and p = 0.75, respectively (with dof of 9416 for both).The best-fit results are given in Table 1.No significant difference is found between the results with p = 0.5 and p = 0.75; all the best-fit values are consistent within the statistical errors.
We confirm that the cutoffpl component is slightly harder in the pre-peak phase than in the post-peak phase.Since the parameters Γ and E fold are tightly correlated with each other, we investigate their statistical errors in the 2-dimensional space.The bottom panel of Figure 2 shows the confidence contours of these parameters for the slim disk (p = 0.5) case.A significant difference is found between the two phases at the 90 % confidence level.A similar result is obtained for the p = 0.75 case as well.We have performed phase-resolved spectral analysis of the broadband X-ray spectra of NGC 5907 ULX1 in its bright phase in 2014.The spectra are well-described by the model consisting of a time-invariant accretion disk component and a pulsed emission component with foreground extinction.The disk component can be equally well described by a geometrically-thin (standard) disk and a geometrically-thick (slim) disk, although the latter would be favored from the nature of the ULX.The pulsed emission component, associated with the corotating accretion flow onto the neutron star within the magnetosphere, can be reasonably modeled by a power law with an exponential cutoff.Its spectral shape, characterized by the combination of the photon index and cutoff energy, is different between the pre-peak and postpeak phases; the former spectrum is harder than the latter one.
The spectra we have presented in this paper were studied in some previous work as well . Fürst et al. (2017) analyzed phase-averaged spectra for the purpose of investigating the spectral changes as a function of the 78day super-orbital phase of the binary system.They first attempted a model similar to ours (i.e., multicolor disk blackbody + cutoff power law), and obtained the best-fit innermost disk temperature consistent with our results (i.e., T in ∼ 0.3 keV).However, this model left systematic residuals above ∼ 10 keV, unlike our results based on the phase-resolved spectroscopy.In their work, this hard excess was attributed to the Comptonization of disk photons, and a simple Comptonization model (simpl: Steiner et al. 2009) was introduced as a convolution component applied to the diskpbb component.As a result, they obtained an innermost disk temperature (T in ) of ∼ 3 keV, significantly higher than our result.We note that the simpl model in XSPEC is an empirical model of Comptonization that converts a fraction of input seed photons into higher-energy photons with a power law distribution.In this model, the spatial distribution of the Comptonizing electrons is assumed to be uniform around the seed photon emitting regions, so all the seed photons are scattered at the same probability (Steiner et al. 2009).This means that the simpl*diskpbb model introduced by Fürst et al. (2017) represents the situation where the accretion disk is uniformly surrounded by the Comptonization region.Given the presence of the pulsation, we believe that the hard X-ray photons originate from the region much closer to the neutron star, and thus this empirical model is not introduced in our analysis.Walton et al. (2018a) conducted phase-resolved spectroscopy, but the phase dependence in the spectral shape was not taken into account.In fact, they extracted spectra from the brightest and faintest quarters of the pulse cycle, and subtracted the latter from the former (i.e., pulse-on -pulse-off) to determine the spectral shape of the pulsed emission component using a cutoffpl model.Subsequently, the phase-resolved spectra were fitted with the diskpbb and cutoffpl components, similar to our analysis, but the photon index and cutoff energy of the pulsed component were fixed to the values derived from the analysis of "pulse-on -pulse-off" spectrum.The innermost disk temperature was then obtained to be ∼ 3 keV, similar to the result of Fürst et al. (2017) based on the simpl * diskpbb model.In our analysis, we have taken into account the spectral shape variation among the phase bins, obtaining a substantially lower T in value (∼ 0.3 keV).Moreover, no systematic residuals are found in the hard X-ray band over the model consisting of the diskpbb plus cutoffpl components.Keeping them in mind, we discuss the implication of our results in the following subsections.

Disk Inclination Angle
As discussed above, the disk emission is well reproduced by either the standard or slim disk model in our spectral analysis.Its bolometric luminosity is obtained to be L standard = (1.9 ± 0.1) × 10 40 erg s −1 for the standard disk case and L slim = (1.7 ± 0.1) × 10 41 erg s −1 for the slim disk case, assuming the distance to the source to be 17.1 Mpc (Tully et al. 2016).Given the extremely high luminosity observed in this source, it is reasonable to assume that the innermost disk is geometrically-thick due to the super-critical accretion, but a geometricallythin part should also be present in the regions outside the spherization radius.Therefore, the actual luminosity of the whole accretion disk is considered to be in between L standard and L slim .On the other hand, the luminosity of the power-law component is obtained to be L pow = 1.0 × 10 41 erg s −1 , regardless of the applied disk model.The uncertainty to this value is negligibly small.The bolometric luminosity of the whole system is, therefore, constrained to be 1.2 ≤ L bol [10 41 erg s −1 ] ≤ 2.8 . (1) The accretion disk is truncated at the magnetospheric radius R M , where the magnetic field pressure of the neutron star becomes comparable to the ram pressure of the accreting matters.Following Equation S6 of Israel et al. (2017b) and assuming that the geometry of the innermost accretion disk is nearly spherical, we estimate this radius to be where B NS , M NS , and R NS are the surface polar magnetic field strength, mass, and radius of the neutron star, respectively.
Note that this equation is based on the analytical model originally developed by Lamb et al. (1973), where a moderate accretion rate is assumed.It is therefore not obvious if this equation is valid for the system with a super-critical accretion rate.Nevertheless, we first estimate the disk inclination angle based on this model for a direct comparison with the previous work.
Given that NGC 5907 ULX1 was in its bright phase at the time of the observations, R M is considered to be smaller than the corotation radius R c , where the Keplerian orbital angular velocity is equal to the angular velocity of the neutron star rotation.Otherwise, the centrifugal force prevents the accretion flow so that the system becomes X-ray dim.Using the spin period of the neutron star P , we derive the corotation radius as: where G is the gravitational constant.Since P = 1.137 s for NGC 5907 ULX1 at the time of the observations, Eq. 3 gives R c = 1.8 × 10 3 km, when M NS = 1.4 M ⊙ is assumed.Therefore, the condition of R M < R c leads to the relationship between L bol and B NS as: This gives a more stringent lower limit on the luminosity than that in Equation 1 when B NS > 6.5 × 10 13 G.If B NS exceeds 9.9×10 13 G, however, there is no luminosity range that satisfies Equations 1 and 4 simultaneously.Therefore, this value is regarded as the upper limit of the magnetic field strength at the neutron star surface.Since the diskpbb model is an approximation of a multi-color accretion disk, its normalization (N ) provides an estimate of the inner-disk radius (R in ) as: (5) (e.g., Makishima et al. 2000), where d 10 , ϵ, κ and θ are the distance to the object in units of 10 kpc, geometric factor, color-hardening factor, and the inclination angle of the disk, respectively.The geometric factor ϵ compensates for the difference in the innermost edge boundary condition between the assumption in the diskpbb model and the actual accretion disk.In the case of the standard disk, ϵ = 0.412 is often applied (Kubota et al. 1998).On the other hand, Vierdayanti et al. (2008) argued that ϵ = 0.353 is more appropriate for the super-critical accretion case, and Soria et al. (2015) indeed used this value to estimate the black hole mass of M83 ULX-1.It is more difficult to estimate an appropriate ϵ value for a magnetized neutron star, because the accretion disk is truncated due to the presence of the magnetosphere.Nonetheless, we assume ϵ = 0.353 in the following calculations.The color-hardening factor κ is introduced to convert the color temperature determined from the spectral fit into the effective temperature of the disk.At sub-Eddington accretion rates, κ ≈ 1.7 is known to be a good approximation (Shimura & Takahara 1995;Davis et al. 2005).However, the value can be as high as 2.0-3.0 at higher accretion rates (e.g., Watarai et al. 2000;Isobe et al. 2012).
As already mentioned, the accretion disk is thought to be truncated at R M , allowing us to assume R in = R M .Therefore, Equations 2 and 5 and the best-fit N value give a constraint on the disk inclination angle as: Figure 3 illustrates the allowed inclination angle (in blue regions) as a function of bolometric luminosity for different values of κ and B NS .The luminosity ranges that are ruled out from Equations 1 and 4 are indicated as the gray regions.We find that a relatively low inclination angle is favored in general, and that a narrower allowable range is obtained for higher κ and/or lower B NS .It should also be noted that, for κ = 3, the mathematically obvious requirement of cos θ ≤ 1 is not satisfied with any combinations of B NS and L bol that are allowed in Equation 4. We conclude, therefore, that moderate values (2.0-2.5) are favored for the color-hardening factor, despite the high accretion rate achieved in this system.Our result is in contrast to the previous study by Fürst et al. (2017), where a higher inclination angle of ∼ 89 • was claimed.The origin of this discrepancy can be explained as follows.In the previous study, the innermost disk temperature T in was obtained to be ∼ 3 keV, significantly higher than our measurement (∼ 0.3 keV).Since the disk luminosity is proportional to T 4 in , a lower normalization value (N ∼ 1 × 10 −3 ; see Table 1 for comparison) was derived to explain the observed X-ray flux.Consequently, the nearly maximum inclination angle (hence cos θ ≈ 0) was obtained from the equation of R M and R in (see Equations 2 and 5).Our result can more naturally explain the detection of the pulsation, because the pulsating accretion flow can hardly be seen in the edge-on view due to the obstruction by the geometrically-thick disk that should be present in the super-critical accretion system.It is also worth noting that no evidence for an eclipse has been detected from this source, which also implies the low-inclination view of the system.
Lastly, we introduce a model developed more recently by Gao & Li (2021) to derive similar constraints on the disk inclination angle.This model was constructed by expanding the magnetic threaded disk models of Ghosh & Lamb (1979a,b) and Wang (1987Wang ( , 1995) ) to the case of super-critical accretion, for the purpose of estimating the magnetic field of several observed ULXPs, including NGC 5907 ULX1.We follow Equation 17 of Gao & Li (2021) to obtain the magnetospheric radius as: where ξ is a factor expressing the magnetospheric radius as a fraction of the Alfvén radius, and L Edd is the Eddington luminosity.A remarkable difference from Equation 2 is that the bolometric luminosity L bol is no longer contained in the equation.Assuming R in = R M (i.e., using Equations 5 and 7), we obtain a constraint on the disk inclination angle as: The model of Gao & Li (2021) also allows us to estimate the magnetic field self-consistently.For NGC 5907 ULX1, they used the spin period and other physical quantities reported in Israel et al. (2017b) (i.e., the 2014 observations that have been revisited in this work) and obtained the possible combinations of B NS and ξ as follows: (B NS [G], ξ) = (6.45× 10 11 , 0.88), (6.41 × 10 11 , 0.88), (1.78 × 10 13 , 0.65), and (3.30 × 10 13 , 0.54); see Table 2 of Gao & Li (2021).The first two combinations indicate low magnetic field, yielding no solution for the inclination angle in Equation 8.The other two combinations with the higher magnetic field allow a range of inclination angles as illustrated in Figure 4. Similar to the previous discussion, a relatively low inclination and moderate κ values are favored.

Accretion Geometry within the Magnetosphere
Within the magnetosphere, the accreting matters follow the magnetic field lines toward the magnetic pole of the neutron star (Basko & Sunyaev 1975, 1976).In this section, we discuss the possible structure of the accretion flow that is responsible for the properties observed in the phase-variable cutoff power-law component, i.e., a harder spectrum observed in the pre-peak phase.
In the case of super-critical accretion, the optically thick region within the magnetosphere is thought to become geometrically broad, presumably spreading across the regions between the neutron star's magnetic pole and the innermost radius of the accretion disk.There are two possible radiation mechanisms that are responsible for the pulsed power-law component: Comptonization of soft photons originating from the neutron star surface and/or the inflowing gas (e.g., Inoue ).In the former case, the photon index and cutoff energy of the power-law component can be considered to represent the average number of scatterings and electron temperature, respectively.In the latter case, on the other hand, these parameters would be related to the apparent temperature gradient and maximum temperature of the accretion flow.In either case, the parameters of the cutoff power law component can be considered to reflect the temperature structure of the accretion flow within the magnetosphere.Therefore, we can attribute the phase-variable spectral properties to an axial-asymmetric temperature distribution of the accretion flow such that a relatively hot region emerges in the pre-peak phase.
Here we propose a possible geometry of the accretion flow as illustrated in Figure 5.We assume that the rotation axes of the neutron star and the accretion disk are well aligned with each other, and that the system is observed in the low inclination angle (e.g., θ ≈ 30 • ).When the magnetic axis is oriented toward the observer, the largest cross section of the magnetosphere is achieved and thus the pulse profile peaks.Since the neutron star is in the spin-up phase ( Ṗ < 0), the magnetic field lines could be dragged in the direction of rotation by the accretion disk (red arrows in Figure 5).As a result, the regions closer to the magnetic pole (where the temperature is relatively high) can be more easily observed in the pre-peak phase than in the post-peak phase, making the spectrum slightly harder in the former.If this

Before pulse peak
After pulse peak Magnetosphere Rotation axis M a g n e t ic a x is Gravitational potential hypothesis is true, similar spectral behavior might be observed in other ULXPs in the spin-up phase.Therefore, a systematic study of ULXPs based on the phaseresolved spectral analysis is highly encouraged.Future numerical work is also urged to confirm if the proposed geometry is indeed realized and the observed spectral properties can be quantitatively reproduced.

CONCLUSIONS
In this paper, we have presented the timing and spectral analysis of the broadband X-ray data of NGC 5907 ULX1 obtained from the simultaneous observations with XMM-Newton and NuSTAR during its ultraluminous phase in July 2014.The spin period and its derivative are determined to be P ∼ 1.1 s and Ṗ ∼ −5×10 −9 s s −1 , respectively.The energy-dependent phase shift is observed in the folded light curves, confirming the previous results of Israel et al. (2017b).
The emission from the accretion disk is reproduced by a multicolor disk blackbody model with an innermost disk temperature of ∼ 0.3 keV.Our measurement indicates a relatively low inclination angle of the binary system, in contrast to the previous claim by Fürst et al. (2017).Our result naturally explains the detection of the pulsation as well as the absence of an eclipse phase in this system.
The spectrum of the pulsed emission component is modeled by a power law with an exponential cutoff.Its photon index and cutoff energy can be considered to reflect the temperature structure of the optically-thick accretion flow inside the magnetosphere.Our phaseresolved spectral analysis indicates that this component is slightly harder in the pre-peak phase than in the postpeak phase.This suggests that the magnetosphere has an asymmetric geometry because of the dragged magnetic field lines, and the regions closer to the magnetic pole emerge in the pre-peak phase.A systematic observations of other ULXPs as well as numerical work are highly encouraged to confirm if this hypothesis is true.APPENDIX A. SPECTRA FROM ALL EIGHT PHASES AND CORRESPONDING BEST-FIT PARAMETERS Figure 6 shows the phase-resolved spectra obtained in the manner described in Section 3.2.In Table 2 we present the complete best-fit parameters.The method of the analysis from which these parameters are obtained is also described in section 3.2.Phase Γ E fold (keV) F (10 −12 erg cm −2 s −1 ) p = 0.5 0.0 − 0.1 0.82 +0.17 Notea Uncertainties were obtained using XSPEC steppar command in the 2D-space.
Figure 1.Pulse profiles from the individual instruments with different energy ranges.For (a), the energy ranges are as follows (from the top to the bottom): 0.2 − 12 keV, 0.2 − 3 keV, 3 − 7 keV and 7 − 12 keV.For (b) and (c), the energy ranges are as follows in the same order: 3 − 30 keV, 3 − 7 keV, 7 − 12 keV and 12 − 30 keV.The profiles are normalized by the average count rate over the observations.

Figure 2 .
Figure 2. (Top) Response unfolded spectra of NGC 5907 ULX1.The pre-peak (phase 0.1 − 0.3) and post-peak (phase 0.6 − 0.8) data are shown in red and blue, respectively.The spectra are rebinned, and those from the different instruments share the same marker for visual clarity.For the same reason, the NuSTAR FPMB data are omitted from the plot.The lines represent the best-fit diskpbb+cutoffpl model components with p = 0.5 (see text).The colored dashed and solid lines correspond to the contribution of the total and cutoffpl models, respectively.The black solid line represents the diskpbb component (for the EPIC-pn only).(Bottom) Confidence contours of the photon index versus the cutoff energy of the cutoffpl model component for the case of p = 0.5.Confidence levels from the pre-peak are shown in red, while those from the post-peak are in blue.The dashed and solid lines represent the 68% and 90% confidence levels, respectively.The circle in each color indicates the best-fit value.
Summary of the Results and Comparison with Preceding Work

Figure 3 .
Figure 3.The inclination angle of the accretion disk as a function of the bolometric luminosity under different assumptions of the color-hardening factor κ ([2.0, 2.2, 2.4]) and magnetic field BNS ([3.0, 5.0, 7.0] × 10 13 G).The blue regions indicate the allowed range of the inclination angle with 1σ confidence level derived from Equation6.The gray areas represent where the luminosity is outside the limits.

Figure 4 .
Figure 4.The inclination angle of the accretion disk as a function of the color-hardening factor κ under different pairs of the magnetic field BNS and ξ derived in Gao & Li (2021).The blue and green regions indicate the allowed range of the inclination angle with a 1σ confidence level derived from Equation 8.

Figure 5 .
Figure 5. Schematic illustration of the possible accretion geometry within the magnetosphere.The red arrows indicate the rotational direction of the accretion disk surrounding the magnetosphere.The gray-scale gradient represents the depth of the gravitational potential.

Figure 6 .
Figure 6.Unfolded spectrum from each phase.XMM-Newton EPIC-pn data are shown in black, NuSTAR FPMA data are shown in red, and FPMB data are shown in blue.The best-fit diskpbb+cutoffpl model with p = 0.5 is superimposed.The dashed lines and the solid lines in each color represent the total model and the cutoffpl component, respectively.The black dotted lines represent the diskpbb component.

Table 1 .
Best-fit Model Parameters for Pre-peak and Post-peak Spectra

Table 2 .
Best-fit Parameters of Pulsed Emission Component for All the Spectra