Decomposing the Spectrum of Ultra-Luminous X-ray Pulsar NGC 300 ULX-1

In Figure 1, we present background subtracted-light curves of ULX-1. Those in the larger panel have a length of the overall observation with a bin width of 700 sec, and one in the top left is a small portion (200 sec) of XMM-Newton pn observation with a bin width of 4 sec. We can confirm the presence of coherent variability with the $\sim 31$-sec pulsation period in the 4 sec bin light curve, whereas nothing significant can be found in the others with the longer-time bins. Thus, the data set contains no apparent X-ray variability but from pulsation, making these data sets ideal to study the pure variability of the pulsating component.

Figure 2 represents how the X-ray intensity of NGC 300 ULX-1 varies in terms of the pulse phase and energy. As described in section 1, we folded the entire data set with the initial period and spin-up rate derived by Walton et al. (2018a), which are 31.7183411308 sec (reference epoch is 57738.65732 MJD) and $-5.563\times 10^{-7}$ sec sec^-1, respectively. The pulse profile is relatively sinusoidal, which is consistent with the previous results (e.g., Carpano et al. 2018), and peaked at $\phi=0.5\--0.6$ throughout $0.3\--25$ keV, suggesting a “single zone” hot spot.

Figure 3 presents an energy dependency of the fractional Root Mean Square (RMS) amplitude, $\sqrt{S^{2}/\bar{x}}$, where $S^{2}$ and $\bar{x}$ are the variance and average of the count rate over the pulsation cycle, respectively. The fractional RMS amplitude becomes larger as the energy increases, reaching $\geq 70\%$ above $3$ keV. This is consistent with the previous results that used the same data sets as the present study (e.g., Carpano et al. 2018; Vasilopoulos et al. 2019). Hence, the spectrum of NGC 300 ULX-1 is expected to be dominated by the pulsating emission component at a higher energy band.

Just like those in ULXPs, AGNs in Seyfert 1 galaxies also tend to exhibit featureless X-ray spectra, allowing multiple physical models to reproduce them without any significant statistical differences. Therefore, adding restrictions to models has been crucial, and spectral variability is one of the helpful tools to resolve this model degeneracy.

Noda et al. (2014) introduced a method called Count-Count Correlation with Positive Offset (C3PO), which utilizes the characteristic X-ray variability of AGN to extract a pair of spectra that composes the original spectrum without relying on any physical models. One is from the component that accounts for the variability, and the other is from the unchanging one. The method is based on a correlation between the count rates of two energy bands. One is from a fixed energy range defined as a reference band, and the other is from a test band that is an arbitrary part of the rest. If we plot the test-band count rate against that of the reference band, the data points ought to exhibit a certain locus whose shape depends on how the spectral shape changes in time. The simplest example is when the variable component changes only its intensity. In this case, the data points will form a straight line on the count rate v.s. count rate plot (CCP) plane, and a product of its slope and the count rate in abscissa represents the actual count rate of the variable component at that energy band. Furthermore, if the locus shows any positive offset at the zero point of the reference band, then the source spectrum is likely to contain a component that does not correlate to the variability with a count rate equivalent to that very offset value.

Noda et al. (2014) found that the CCPs at the bright state of a highly variable AGN, NGC 3227, form straight loci, and each of them shows an apparent positive offset if they extrapolate the data down to the count rate where that of the reference band reaches zero. It clearly suggests, as described above, that the X-ray spectral continuum of NGC 3227 consists of at least two components. One is something uncorrelated (or stable) against the variability representing the positive offset in CCP, and the other is the variable component that changes only its intensity in time. The authors fitted each CCP with a linear function expressed as

where $x$, $y$, $a$, and $b$ are the reference-band count rate, test-band count rate, slope, and y-intercept value, respectively. Here, $a$ and $b$ are set as free parameters, and the derived value $ax_{0}$ ($x_{0}$ is the average of $x$ used in the fitting) and $b$ represent the actual count rate of the variable and stable component at that energy band, respectively. By repeating this procedure in the other energy bands, Noda et al. (2014) successfully extracted the spectrum of the stable component with a thermal emission shape and a variable one with an absorbed power-law figure hidden under the featureless X-ray spectrum of NGC 3227. In the present study, we apply the same C3PO method as Noda et al. (2014) to NGC 300 ULX-1, but in a pulse phase-dependent manner.

Figure 4 presents the extracted CCP of ULX-1. Taking the high photon statistics of XMM-Newton and the highest pulse fraction (or variability) at $>7$ keV (Figure 3) into account, we select $7.0\--10.0$ keV of XMM-Newton pn as the reference band (horizontal axis). To grasp the overall behavior of the CCP, we maximized the statistics by designating the rest of the entire energy band, $0.3\--7.0$ keV, as the test band (vertical axis). Instead of the raw count rates utilized in Noda et al. (2014), each data point in this CCP represents an average count rate within a divided portion ($1/20=0.05$ cycle per data point in this case) of the $\sim 31.7$ sec pulsation cycle. In short, the CCPs in the present paper reflect how the count rate in each energy band varies as a function of the pulsation phase.

The CCP of ULX-1 is forming a rather straight locus with a “break” at $\sim 0.03$ count sec^-1, which corresponds to a pulse phase of $\phi\sim 0.4$ and $\phi\sim 0.75$. It suggests that the variable component, emission bound to the pulsating accretion flow, increases its X-ray intensity without changing its spectral shape within phase intervals below or above this break. Since the data above the bent are forming a locus with a shallower slope, the appearance of the spectra should be different between these two phases.

A similar breaking CCP is also reported in NGC 3227 by Noda et al. (2014). Instead of explaining the entire data points with a single linear line as Noda et al. (2011) and Noda et al. (2013) did, the authors tested two alternative functions that could fit such curving CCPs. One is a piecewise-segmented linear function that consists of a pair of linear functions individually explaining the data points below and above the break separately. The other is a single power-law function expressed as $y=Mx^{N}$ (e.g., Uttley & McHardy 2005), where $M$ and $N$ are left free to vary. Following Noda et al. (2014), we tested these two possible solutions as shown in Figure 4. The piecewise-segmented linear function is defined around the breaking point $x_{\rm b}$ as,

where $a_{i}$ and $b_{i}$ ($i=1,2$) are the slopes and intercepts of individual linear functions, respectively. Hence, the best-fit slopes and intercepts automatically derive $x_{b}$. In the regression analysis, we utilized the ROOT analysis package developed by CERN. The errors in $x$ are projected to the $y$-axis direction to take contributions from both $x$ and $y$ into account by calculating chi-square at each data point as $\chi^{2}=(y-f(x))/(\sigma_{y}^{2}+\sigma_{x}^{2}f^{\prime}(x)^{2})$, where $\sigma_{x}$, $\sigma_{y}$, $f(x)$, and $f^{\prime}(x)$ are the errors in $x$ and $y$, the fitting function, and the derivative of the function, respectively.

Although the piecewise linear function gave a slightly better fit ($\chi^{2}/$degree of freedom $=29.2/16$) than the single power-law function ($\chi^{2}/$degree of freedom $=34.2/18$), the difference is rather marginal. If we calculate the Bayesian Information Criteria, which can be expressed using the likelihood function $L$, number of parameter $k$, and number of data points $n$ as $-2\ln L+k\ln n$, the difference between the models is $<1$. According to criteria by Kass & Raftery (1995), this is insufficient to claim that the piecewise linear function is statistically preferable over the power-law function.

Unlike Noda et al. (2014), we were unsuccessful in statistically ruling out the power-law function due to the limited number of data points in the CCP. Since the difference between the two functions becomes significant around the breaking point, it may be possible to distinguish one from another by observing ULX-1 with a long exposure or instruments with larger effective areas. We leave this as a future work and adopt the results on piecewise linear function as a working hypothesis of the following analysis. The best-fit value of the breaking count rate is $0.028$ count sec^-1, and hence we hereafter refer to the phase below this count rate ($0.0\leq\phi\leq 0.4$ and $0.75\leq\phi\leq 1.0$) as “faint phase” and above it ($0.4<\phi<0.5$) as “bright phase”.

To extract the spectra of stable and variable components, we divided the $0.3\--7.0$ keV of XMM-Newton pn and $3.0\--25.0$ keV of NuSTAR FPM into 15 and 5 sub-bands, respectively, and created CCPs for each using $7.0\--10.0$ keV count rate of XMM-Newton pn as a reference (Figures 5, 6 and 15). The sub-bands are defined as those consisting of data points that contain at least 20 counts and have logarithmically-equal width within the individual instruments except for the $10\--25$ keV bands of NuSTAR FPMs. Thus, the same breaking feature as Figure 4 is present in the other energy band. The change in slope is apparent in the lower energy bands and becomes more ambiguous as the energy of the test band reaches 10 keV or higher.

As Noda et al. (2014), we fit the data points with count rates below and above the break in each CCP separately with the respective linear functions shown in blue and red solid lines in Figures 5, 6, and 15. The goodness of fit and obtained parameters are summarized in Table 1. The pairs of linear functions successfully reproduced the individual CCPs, and some exhibited non-zero y-intercept values. Hence, following the procedure described above, we can generate two pairs of spectra from the best-fit values for each phase. Combining XMM-Newton and NuSTAR, we obtained a stable-component spectrum with 7 bins in $0.30\--1.61$ keV via data points in the faint phase area and another with 20 bins in $0.3\--10$ keV from those on the other side. As for the variable component, we obtained a spectrum with 25 bins between $0.3\--25$ keV from the faint phase, whereas those from the bright phase produced one with 22 valid bins ranging from $0.56$ keV to $25$ keV.

In a typical phase-resolved spectral analysis, many compare spectra extracted from limited time intervals around the peak and bottom of the pulsation. For example, Koliopanos, F. et al. (2019), who analyzed the same data set as the present study, utilized 60% (30% each for on and off-pulse spectra) of a rotation cycle and left the rest from the phase-resolved analysis. Although this may clarify the difference between the on-pulse and off-pulse spectra, one can also discard the photon statistics excessively by restricting the time interval, and lead several models to degenerate.

In contrast to the ordinary phase-resolved analysis, our study enables us to add further information to untangle the model degeneracy and utilize the data more efficiently. As the CCPs clarified in the previous section, the entire pulsation cycle of NGC 300 ULX-1 can be categorized into two phase intervals, namely the “faint phase” and the “bright phase”. Since we now visually know from the CCPs that the spectral shape is constant within each phase, we do not have to limit the time interval further to see the change in spectrum. Hence, we hereafter divide the entire observational data into these two phases and extract spectra from each using all of the data therein to study what spectral model composition can explain the spectrum from each epoch with as high photon statistics as possible.

We present the bright phase (black) and the faint phase (red) spectra in the top panel of Figure 7. To tentatively unfold the spectra with the instrumental response, they are shown in ratios over a power-law model with a photon index of $\Gamma=2$. Thus, ULX-1 exhibits a hard ($\Gamma\sim 1.4$) and continuum-dominated spectrum throughout the pulsation phase. As the pulsating component decreases its intensity, an additional thermal component peaking at 1 keV becomes more apparent above the continuum. While the ratio between the two spectra (the bottom panel of Figure 7) indicates a spectral hardening in the $0.3\--5$ keV band, not as apparent in the $>5$ keV band. This is consistent with the behavior we saw in the “breaking” feature of the CCPs presented in section 3.2. In the following section, we study the underlying emission compositions of these two spectra along with those newly extracted via the C3PO method.

In this section, we give a possible implication to the origin of the convex non-pulsating spectral component extracted from the faint phase. According to theoretical studies on super-critically accreting objects, the entire accretion flow, from mass donating star to accretor, can be roughly separated into three regions in terms of distance from the center. One is the most outer region of the accretion flow where the radiative cooling is efficient enough to form an optically thick and geometrically thin standard accretion disk. In this region, the disk emits a multi-color disk blackbody spectrum with a radial surface temperature dependence of $\propto r^{-0.75}$ (Shakura & Sunyaev, 1973).

The closer to the center the stronger the emission becomes, and at a certain radius, the radiation pressure eventually overwhelms the self-gravitational pull of the disk. This forms the second flow region, in which the accretion disk starts to puff up and deviate from the ordinary standard disk regime. Within this radius, the accretion flow is now advection-dominated, in which a significant fraction of generated photons inside the disk are engulfed instead of being radiated from the disk surface. Although the region also emits a multi-temperature blackbody as the outer “standard disk” region, the radial temperature dependence is expected to be flatter as briefly described in section 3.3, namely $\propto r^{-0.5}$ (e.g., Watarai et al. 2000, 2001). In addition, a cool-optically-thick wind is expected to be launched from this region due to the intense radiation pressure (e.g., Ohsuga et al. 2003; Kawashima et al. 2012). The edge of this outflow can form a photosphere that emits an additional optically thick emission with a characteristic temperature.

If the central object is a strongly magnetized neutron star, its magnetic pressure eventually overcomes the gas pressure of the flow at a certain point closer to the central object. Hence within this radius, the magnetic force restricts the accreting matters to move only along the magnetic field. This is the third characteristic region in the accretion flow, and it is called an accretion column (e.g., Basko & Sunyaev 1975, 1976). As the neutron star rotates, the entire column precesses along the magnetic field, accounting for the pulsating emission component.

Since the stable component is unrelated to the pulsation, the emission clearly originates from some regions outside the accretion column, where the dipole magnetic field of the neutron star is too weak to capture the infalling matters. Therefore, it should be emission from either the outer or inner region of the accretion disk. According to the spectral analysis in section 3.3.2, two types of thermal emission models, diskbb and diskpbb, gave acceptable fits. Since the single-temperature blackbody model failed to reproduce the spectra, we consider it unlikely that a photosphere of the outflow is the origin of the stable component. Although allowing the radial temperature profile to vary (using diskpbb) did not improve the fit and we could not statistically distinguish whether the data favor the disk to be in the standard or advection-dominated regime, we hereafter assume the former as a working hypothesis and proceed to further discussion.

Since diskbb is an approximation of the standard accretion disk, its normalization gives an apparent inner-disk radius. The realistic radius $R_{\rm in}$ can be calculated by using the model normalization parameter $N$ as

(e.g., Makishima et al. 2000). Here, $D_{10}$, $\xi$, $\kappa$, and $\theta_{\rm i}$ represent the distance to the object in units of 10 kpc, a correction factor, color-hardening factor, and the inclination angle of the disk, respectively.

The factor $\xi$ corrects for the difference in the inner-most edge boundary condition between diskbb and the standard accretion disk. The value is known to be $\xi=0.412$ (Kubota et al., 1998) for a realistic standard accretion disk model that takes general relativity effects around a black hole into account. Due to the presence of the interrupting magnetic field, it is rather complex and challenging to estimate $\xi$ for the disk around an accreting neutron star. Some theoretical works suggest that the accretion disk may resemble the $\xi=0.412$ case, still, the value strongly depends on the details of the accretion flow (e.g., Nixon & Pringle 2021). Since $\xi$ for accreting neutron stars is thus unclear, here we assume the most popular value, $\xi=0.412$ (Kubota et al., 1998), as other studies on accreting neutron star low mass X-ray binaries (e.g., Sakurai et al. 2014)

The color-hardening factor $\kappa$ is a ratio of the color temperature to the effective temperature. The value is known to be in between $\kappa=1.5\--2.0$ (Shimura & Takahara, 1995) at the sub-Eddington regime and recent numerical study suggests that the factor has a weak dependency on the mass accretion rate up to Eddington (e.g., Davis & El-Abd 2019). Although the total mass accretion rate is in the a super-Eddington regime for a neutron star, here we assume that the local accretion mass density is still in the sub-Eddington level for the accretion flow at the radius where the disk is formed. Hence, we employ the most commonly used value of $\kappa=1.7$ (Shimura & Takahara, 1995).

If we assume $D$ to be the distance to NGC 300 (1.9 Mpc; Gieren et al. 2005), then the present results in Table 2 and Equation 4 give a radius of $R_{\rm in}=720^{+220}_{-120}/\sqrt{\cos\theta_{\rm i}}$ km. Although the inclination angle of this system is still unknown, considering the detection of pulsation and the fact that no evidence of eclipse is present, we may assume that the system is nearly face-on to the observer ($\theta_{\rm i}\sim 0$). As discussed in section 4.1, $R_{\rm in}$ can imply the radius where the accretion flow shifts its behavior from the standard accretion disk to an accretion column; the magnetospheric radius $R_{\rm M}$.

In addition to the inner-disk radius, we can calculate the mass accretion rate at the radius by utilizing the theory of standard accretion disk (Shakura & Sunyaev, 1973) as

where $G$, $M$, and $L_{\rm disk}$ are the gravitational constant, central object mass, and disk luminosity, respectively. Assuming a typical neutron star mass of $M=1.4M_{\rm\odot}$ and substituting values for $L_{\rm disk}$ and $R_{\rm in}$ from the present result, equation 5 gives a mass accretion rate of $\dot{M}\sim 1.3\times 10^{20}$ g sec^-1, which is $\sim 50$ times the Eddington rate of $1.4M_{\rm\odot}$ neutron star ($\sim 2.8\times 10^{18}$ g sec^-1). In the following sections, we derive several physical values from these $R_{\rm in}$ and $\dot{M}$ and give a plausible explanation for the observed characteristics of NGC 300 ULX-1.

Utilizing the obtained observational values, we estimate the strength of the dipole magnetic field of ULX-1 by following similar discussion as Walton et al. (2018b). The magnetic field of a mass-accreting neutron star with a certain mass accretion rate $\dot{M}$ can be expressed as

(Ghosh & Lamb, 1979b; Lai, 2014; Fürst, F. et al., 2017). Since $R_{\rm in}$ is now observable and $\dot{M}$ can be calculated via equation 5 utilizing observed $L_{\rm disk}$, we are able to derive the magnetic field $B$ by assuming a typical neutron star mass of $M=1.4M_{\rm\odot}$ from equation 10.

In Figure 13, we present the relation between the mass accretion rate and the magnetic field derived from equations 10 (black) and 5 (red). Each curve is drawn using the result obtained in the spectral analysis. According to this figure, our estimation of the magnetic field strength of ULX-1 lies in $2\--7\times 10^{12}$ G (from the lowest to the highest tip of the red-hatched region).

The estimated value is consistent with those made by other methods. For example, Vasilopoulos et al. (2018) estimated the strength of field as $5\times 10^{12}$ G by utilizing the observed spin period evolution and a theory based on induced torque to neutron stars from accreting matters. Walton et al. (2018a) have found a possible absorption line feature that may be interpreted as a cyclotron resonance scattering feature. Assuming the feature originates from electron scattering, the authors concluded that the estimated magnetic field is $\sim 10^{12}$ G.

Finally, let us discuss the possible structure of the pulsating accretion flow in NGC 300 ULX-1 from the variability we saw in the present analysis. As described in section 3.2 and section 3.3, we have revealed that the variability of the pulsating component of this ULXP can be divided into at least two phases, the faint phase ($0.0<\phi<0.2$ and $0.5<\phi<1.0$) and the bright phase ($0.2\leq\phi\leq 0.5$).

In the faint phase, the pulsating component exhibited a hard continuum ranging from 0.5 keV to 25 keV with a rollover at $\sim 7$ keV. Although we could roughly explain the spectrum with models that emit photons in a wide energy band with a characteristic cutoff temperature (multi-color disk blackbody, cutoff power law, and Comptonization), it required an extra component to explain the extending continuum at $>10$ keV. The result is consistent with the previous studies of ULXs including the same ULXP (e.g., Carpano et al. 2018; Walton et al. 2018b; Koliopanos, F. et al. 2019), and we employed an extra model to resolve this discrepancy. Simply adding a power-law model could not solve the residual at $>10$ keV because the model extends to infinity in both energy directions and causes a conflict with the stable component that dominates the flux in the lower energy band. Thus, we concluded that the data favor a model that sharply drops off at the lower energy (e.g., Rayleigh-Jeans of the blackbody), and a model consisting of Comptonization plus its power-law scattering fits the data the best so far.

As discussed in the previous section, strongly magnetized neutron stars can form a magnetosphere, which forces the accreting matter to fall along its magnetic field. The magnetosphere creates a cylindrical-shaped accretion flow at a close region to the magnetic pole so-called accretion column (e.g., Basko & Sunyaev 1975, 1976). Within this column, free-falling matters are shock heated at a particular height from the stellar surface and emit high-energy X-ray photons. If the accretion rate gets close to the Eddington limit, these generated photons begin to escape from the sidewalls of the column rather than the direction of the magnetic field, creating an emission pattern perpendicular to the magnetic field called a “fan beam” (Basko & Sunyaev, 1976). Some of these fan-beam photons can irradiate an optically thick region trapped around the boundary of the magnetosphere. The incident photons will experience photoelectric absorption and multiple Compton scatterings within this region, and some eventually escape the system as reprocessed thermal emission (e.g., Mushtukov et al. 2017). Other fractions of the fan beam photons may instead irradiate the neutron star surface close to the magnetic pole and be reflected as a secondary emission pattern called a “polar beam” (e.g., Trümper et al. 2013; Poutanen et al. 2013). In this case, photons are re-emitted parallel to the magnetic field with a spectrum harder than the original fan beam. Thus, the accretion flows around highly magnetized mass accreting neutron stars can form multiple emission regions, and several observations on other ULXPs (e.g., Koliopanos & Vasilopoulos 2018 for SMC X-3, and Bykov et al. 2022 for Swift J0243.6+6124) have supported such a complicated emission geometry. Although the sources were not in the ultra-luminous state, the resent X-ray polarization observations on accreting neutron star binaries (Tsygankov et al. 2022 for Cen X-3, and Marshall et al. 2022 for 4U 1626-67) have also hinted the presence of the multi-zone pulsed emission.

Some numerical studies suggest (e.g., Mushtukov et al. 2017, 2019) that at well above the Eddington rate, the reprocessing region can extend further down to the neutron star to obscure the central hard X-ray emitting region. As almost the entire region within the magnetosphere has shifted to the reprocessing area, the accretion flow at this rate is no longer regarded as a column but rather as a “curtain” (for example, see Figure 1 in Mushtukov et al. 2017). It shields and reprocesses most of the hard X-ray photons from the central fan beam into thermal black-body-like emission. Since this reprocessing optically thick curtain has a particular temperature gradient along the magnetic field, one can approximate its overall spectrum with the summation of black bodies from respective temperature regions. A numerical study estimates the lowest and highest temperature of the black bodies to be around sub keV and a few keV for a neutron star accreting at $5\times 10^{39}$ erg sec^-1 with a magnetic field of $10^{13}$ G (Mushtukov et al., 2017). The estimated value is roughly consistent with the present observational result. The variable component spectrum exhibited breaks at both ends in energy, and we successfully explained the lower one with the Rayleigh-Jeans of 0.16 keV blackbody. Although the Comptonization model, not the multi-color blackbody as expected in the theoretical study, gave the best fit to the data, we consider that the model just happened to be the one to approximate the actual temperature gradient of this accretion flow. The smooth single-peaked pulse profile (Figure 3) suggests that the emission from the other side of the pole is not reaching the observer, which also supports the idea that the reprocessing region covers most of the magnetosphere.

According to several theoretical studies, including numerical simulations (e.g., Ohsuga et al. 2005; Ohsuga 2007; Ohsuga & Mineshige 2011; Kawashima et al. 2012; Jiang et al. 2014; Sadowski et al. 2014), an extreme accretion rate as ULXPs may derive the surrounding disk to launch optically thick and cool outflows. Since this creates a tall funnel-like structure around the central source, it can easily interfere with the pulsating emission and change its spectral shape by Compton down scattering or absorbing photons as the entire accretion curtain precesses around the spin axis. Kosec et al. (2018) reported the presence of blue-shifted absorption lines in the same data set with a $\sim 3\sigma$ significance. The detection is evidence of highly-ionized matters outflowing with $\sim 20\%$ of the speed of light. However, the spectrum of the pulsating component does not change its shape as the intensity varies within each phase interval (see the CCPs in section 3.2). The behavior indicates that the emission appearance of the individual accretion curtain regions, which accounts for the intensity change in each phase interval, is somehow uniform and does not depend on the rotation phase; suggesting only the solid angle of the emission region to the observer changes as the curtain rotates with the neutron star. Furthermore, (Kosec et al., 2018) estimated the column density of the outflow to be $1.2^{+1.9}_{-0.6}\times 10^{23}$ cm^-2, which is relatively transparent to the Thomson scattering ($<10^{24}$ cm^-2). According to several numerical simulations of super-critical accretions, such Thomson-thin gas can leak out from the accretion curtain (Abolmasov & Lipunova, 2022) and be accelerated up to $10\--40\%$ of the speed of light via intense radiation pressure and magnetic reconnections near the magnetic pole (Takahashi & Ohsuga, 2017). Considering the observational and simulation results above, we propose that, instead of the disk, the intense radiation of the rotating accretion curtain is blowing the thin matter leaked out near the neutron star away. Since the strong magnetic field is likely truncating the disk before it derives the geometrically/optically thick super-critical winds (see section 4.1.1), the central curtain emission is visible without being interfered with by the matters within the disk.

Although the theoretical studies thus qualitatively describe most of the spectral behavior of the pulsating component below 10 keV, we still do not have a good explanation for the tail component above that. Due to the limited signal-to-noise ratio above 25 keV, it is rather challenging to determine whether the feature above 10 keV originates from an absorption line argued by Walton et al. (2018b) or an extending power-law continuum. We expect a future hard X-ray mission to resolve this problem.

As the source reaches the bright phase (or the pulse peak), the component that used to be variable in the faint phase (blue solid line in Figure 14) halts to increase its intensity, and an alternative variable component (magenta solid line in Figure 14) with a similar but slightly harder continuum emerges. It strongly indicates that the accretion flow accounting for the pulsation consists of at least two representative emission regions. Furthermore, considering that the harder component is observable only in the bright phase ($\sim 30\%$ of the pulsation period), the emission can be slightly collimated toward the magnetic pole axis. Such a multi-region structure has never been reported in the rotating accretion component of ULXPs.

In the top half of Figure 14, we present a schematic cross-section drawing of a possible accretion geometry that explains our observational result. As described above, an optically thick curtain is likely formed along the magnetic field due to the super-Eddington accretion. Since the magnetic field confines the accretion flow into a small region near the magnetic pole, the entire curtain creates a funnel structure, which obscures most of the neutron star and the central hard X-ray (fan beam) emitting region from the observer. Such curtains around magnetized neutron stars are also suggested by theoretical studies including numerical simulations (e.g., Takahashi & Ohsuga 2017; Mushtukov et al. 2019; Abarca et al. 2021; Inoue et al. 2023). The entire curtain exhibits a reprocessed multi-color blackbody spectrum shown in blue lines. The area near the stellar surface (the magenta region) is observable only when the opening cone faces the direction toward the observer as the flow precesses. Thus, the hard variable component shown in the magenta line emerges in a limited phase ($0.2\leq\phi\leq 0.5$), whereas the blue reprocessed component is present in the spectrum at all times.

Since this emerging continuum is harder than that in the faint phase, it can be interpreted as an emission from the polar beam. This is because the polar beam is a scattered-off fan beam component, in which the lower energy part of the original emission experiences a photoelectric absorption at the neutron star surface. The emission of the polar beam typically extends up to $>10$ keV (e.g., Koliopanos, F. et al. 2019; Trümper et al. 2013), whereas that of NGC 300 ULX-1 rolls over at lower energy as $\sim 7$ keV. Since the spectral shape is relatively similar to the variable component in the faint phase, the polar beam can be also reprocessed by the optically thick curtain, and we might have observed the remnant of its original hard spectrum. Recent two-dimensional simulation indicates (Kawashima & Ohsuga, 2020) that the super-critical accretion flow onto a highly magnetized neutron star can form a complex structure inside its opening funnel. Also, pulse shape and pulse fraction depends on the the structure of the funnel and observer’s viewing angle (Inoue et al., 2020). Although the expected spectrum from such a complex structure is yet to be available, the present result might support such a simulation.

Finally, we discuss how the present accretion scenario is related to a longer time-scale spectral evolution. According to the recent SWIFT/XRT and NICER monitoring campaign on ULX-1, the source has continuously decreased its X-ray flux by an order of magnitude since this observation in 2016 (e.g., Ray et al. 2019; Vasilopoulos et al. 2019; Ng et al. 2022). Although the cause of the dimming is still unclear, Vasilopoulos et al. (2019) suggested that obscuration by radiation pressure-dominated disk and its outflow can explain the phenomenon rather than a decrease in intrinsic mass accretion rate. This is because the source kept spinning up at the same rate as that in the brightest era, meaning the central neutron star continuously gained an equivalent amount of accreting matters, i.e., spin-up torque, in the dimming phase. Ng et al. (2022) also supported this scenario by reporting spectral softening and occasional detections of possible partially-covered disk emissions during this X-ray flux decrease. As we calculated in section 4.1.1, the estimated $R_{\rm M}$ of ULX-1 is comparable to or marginally larger than $R_{\rm sp}$ (or $R_{\rm trap}$) in this particular observation. Therefore, a slight increase in mass accretion rate may result in a formation of a geometrically thick radiation pressure-dominated disk region that launches thick outflows. Given that the intrinsic mass accretion rate has been comparable to or possibly higher than the present data, we suggest that a radiation-pressure-dominated disk has formed during the dimming phase due to a further decrease of $R_{\rm M}$ (or an increase of $R_{\rm sp}$) from the 2016 observation and such a thick disk has obscured the central neutron star region as it precesses along the line of sight due to a physical mechanism; such as Lense-Thirring precession (e.g., Middleton et al. 2017) as Vasilopoulos et al. (2019) and Ng et al. (2022) proposed.