NICER and Fermi GBM Observations of the First Galactic Ultraluminous X-Ray Pulsar Swift J0243.6+6124

NICER was launched and installed as an external payload on the International Space Station in 2017 June. NICER contains one instrument, the X-ray Timing Instrument (XTI), sensitive in the 0.2–12 keV range, described in detail in Gendreau et al. (2016) and briefly summarized here. The XTI is composed of an array of 56 X-ray "concentrator" optics, each with an associated focal plane module (FPM) containing a single-pixel detector. Each of the 56 X-ray-concentrator modules consists of 24 nested grazing-incidence gold-coated aluminum foil mirrors, parabolically shaped with a common focal length. Details of the XTI detector subsystem are described in detail in Prigozhin et al. (2016) and summarized here. The FPMs each consist of a silicon drift detector integrated with a custom pre-amplifier. The FPMs, 52 of which are operating on orbit, are arranged into 7 groups of 8. Each group of FPMs is controlled by a set of electronics called a Measurement/Power Unit (MPU) slice. Each MPU slice provides the readout and supporting circuitry for its group of 8 FPMs and is completely independent of the other 6 slices.

The GBM instrument comprises 14 detectors: 12 Sodium Iodide (NaI) detectors sensitive from 8 keV to 1 MeV pointed at various directions in order to cover the whole unocculted sky, and 2 Bismuth Germanate (BGO) detectors sensitive from 200 keV to 40 MeV on either side of the Fermi spacecraft (Meegan et al. 2009). Only the NaI detectors are used for this analysis. GBM data are available continuously, when the spacecraft is outside of the South Atlantic Anomaly (SAA). There are three continuous data types: CTIME (0.256 s, 8 energy channels), used for the public pulse frequency and pulsed flux measurements (shown in Figure 2, with updates found at https://gammaray.nsstc.nasa.gov/gbm/science/pulsars.html) and total flux measurements (https://gammaray.nsstc.nasa.gov/gbm/science/earth_occ.html); CSPEC (4.096 s, 128 energy channels); and CTTE (10 μs, 128 channels) event-by-event data, used for this analysis.

The Swift BAT is a hard X-ray monitor (Krimm et al. 2013) composed of an array of CdZnTe detectors that are sensitive in the 15–150 keV range, with a total detecting area of 5200 cm². It is a coded mask instrument with a 1.4 steradian field-of-view. We use the BAT transient monitor results¹⁸ (15–50 keV), provided by the BAT team, to model the accretion torques applied to the neutron star.

The data analysis procedure for this bright source differs from that used for fainter sources in order to correctly account for detector deadtime. Analysis was done beginning with the unfiltered event files for each MPU slice, and the individual slice data sets were kept separate. HEASOFT version 6.22.1 and NICERDAS version 2018 February 22_V002d were used for this analysis. Good time intervals (GTI) were selected using the FTOOL nimaketime with the following screening criteria: ISS outside the NICER-specific SAA boundary, NICER in tracking mode with pointing direction <0015 from the source direction with at least 38 detectors enabled, source elevation >30° above the Earth limb, and source direction at least 40° from the bright Earth. These GTIs were applied to the unfiltered and uncalibrated events for each MPU slice separately using niextract-events.

Two types of deadtime must be separately accounted for in our analysis. First, event-by-event deadtime occurs because processing each event in each detector takes a certain amount of time and that detector is inactive, hence "dead," while processing an individual event. Second, deadtime can occur because the buffer in each MPU slice reaches its maximum and events are dropped because there is no further room to record them. The first type of deadtime is tracked as a DEADTIME column in the NICER event files and the second (called GTI exposure in this paper) is tracked in the GTI extensions of the NICER event files. Below we account for both types of deadtime in the generation of our NICER pulse profiles and light curves.

The standard NICER calibration was applied to the events from each MPU slice separately, using nimpucal. This tool produces a Pulse Invariant (PI) column by applying the gain correction specified in the calibration database, nixtiflightpi20170601v001.fits, to the raw pulse height measurements. The gain scale is set so that PI = E/10 eV, where E is the nominal photon energy. Both the calibrated and uncalibrated files are barycentered using barycorr with the DE430 JPL ephemeris. The nimpucal tool automatically removes events that cannot be calibrated. Because deadtime arises from all events, processing from this point proceeds in parallel for each MPU slice, with the calibrated files for each MPU slice being used to accumulate X-ray events and GTI exposure in the pulse profiles while the uncalibrated files are used to accumulate deadtime.

Pulse profiles were generated for each observation using the Fermi GBM orbital ephemeris described in Section 4.1, along with the spin frequency and frequency derivative measured with GBM on the same day as the NICER observation, and corrected to the mid-point time of the NICER observation. This process was done in a series of steps described in the following paragraphs.

First, source events were selected to form the pulse profiles. Starting with the barycentered, calibrated, but unfiltered event files for each MPU slice, events were selected that were in the chosen energy range (e.g., PI = 20–1200 for 0.2–12 keV for the full profiles; PI corresponding to 0.2–1, 1–2, 2–3, 3–5, 5–8, and 8–12 keV for the individual energy bands; and PI corresponding to 4–7 and 7–10 keV for the hardness ratio comparisons with Reig & Nespoli 2013). Forced triggers, overshoot, and undershoot events (those events with bits 5, 6, or 7, respectively, set in the EVENT_FLAGS column) were excluded. After correcting the event times for the pulsar's orbit, the pulse phase was computed and events were accumulated in each of 100 phase bins across the pulse period. These folded event profiles were accumulated separately for each MPU slice.

Second, the GTI exposure for each individual phase bin was then computed, again separately for each MPU slice. These GTIs were taken from the GTI extension of the calibrated event files. In the brightest observations, shown in Figure 3, the exposure time per phase bin varies dramatically from bin-to-bin. This is caused primarily by buffers filling up at each MPU slice, which results in very short GTIs. To account for these bin-to-bin exposure time variations, the GTIs longer than a phase bin were divided into pieces 1/10 of a phase bin long. Their start times were corrected for the pulsar's orbit, a pulse phase was computed for each time, and the total exposure was accumulated for each phase bin. GTIs shorter than a phase bin were corrected and accumulated without slicing.

Zoom In Zoom Out Reset image size

Third, deadtimes were obtained from the unfiltered event files in order to include all events that could contribute to deadtime, both the good events included in the pulse profile and the non-source events filtered out of the profile. Deadtime for each event is stored in the NICER fits files as a value from 0 to 127. Each of these values V_dead corresponds to a range of 16 clock ticks (i.e., 0 is 0–15 ticks, 1 is 16–31, and so on), with a maximum recorded deadtime of 127 or 2032 ticks. The deadtime is computed in seconds as t_dead = (V_dead + 0.5) × 16/(25.8 × 10⁶ Hz), where the factor of 0.5 shifts to the center of the recorded interval, 16 is the number of ticks per interval, and 25.8 MHz is the clock frequency. For each event, the event times were corrected for the pulsar's orbit, phases were calculated, and the deadtime was accumulated in each phase bin. Like the folded events and folded GTI intervals, the deadtimes were accumulated for each MPU slice separately.

Fourth, to obtain the deadtime corrected exposure time for each phase bin, , the total deadtime for each phase bin was subtracted from the GTI exposure for each bin. This resulted in an exposure time for each of 100 phase bins for each MPU slice. Table 1 lists the total deadtime corrected exposure averaged over the seven MPU slices for each observation. The exposure is typically the same for all MPU slices to within 1% for average count rates <1000 counts s⁻¹. For the highest count rates, the exposure can vary by as much as 10% from MPU slice to MPU slice.

Finally, the count rate in each phase bin was obtained by simply dividing the total number of filtered events in each phase bin by the deadtime corrected exposure in each phase bin. Poisson-distributed uncertainties were assumed, and were calculated as the square root of the counts divided by the deadtime corrected exposure in each bin. The last step was to combine the count rate pulse profiles for each MPU slice into a single NICER pulse profile for each observation. Because 52 of 56 NICER detectors are operating on orbit, the MPU slices have differing numbers of inputs. The count rate profiles were corrected by a simple factor so that all would have the expected rate for eight detectors. The profiles from MPU slices 1 and 6 were multiplied by 8/7, and the profile from MPU slice 2 was multiplied by 8/6; then the profiles from all seven MPU slices were added together and their errors were combined in quadrature.

NICER observations of Swift J0243.6+6124 began early in the outburst on 2017 October 2, triggered by an alert from MAXI (Sugita et al. 2017b) on 2017 September 29. The outburst profile is asymmetric, with a faster rise than fall as shown in Figure 1, as is common for giant outbursts in Be/X-ray binaries (Reig & Nespoli 2013). Inset plots in Figure 1 show pulse profile evolution in the full 0.2–12 keV band throughout the outburst. The pulse shapes are similar at similar count rates both in the rise and fall of the outburst, evolving from highly complex (Figure 1 insets (a) and (h)) at lower count rates, to single peaked (inset (b) and (g)), to double peaked and nearly symmetric (inset (c) and (f)), and to double peaked but asymmetric at the outburst peak (inset (d) and (e)). Detailed energy dependent pulse profiles for each NICER observation are shown in the Appendix.

NICER observations continued as the outburst faded. The flux reached a minimum in the NICER band on MJD 58204 (see Table 1). Pulsations were still detected, and the 0.2–12 keV profile was a simple symmetric single peak, similar to MJD 58176, the last profile shown in Figure 8 in the Appendix. Then the flux again began to increase and a normal (type I) outburst was detected with NICER, GBM, and Swift BAT peaking on MJD 58220. A second normal outburst peaking near MJD 58247 was detected with GBM and Swift BAT, but was not observed with NICER.

The root mean squared pulsed fraction f_rms was computed for the full 0.2–12 keV band and for each of six energy bands: 0.2–1, 1–2, 2–3, 3–5, 5–8, and 8–12 keV. The pulsed fraction is given by

$$\begin{eqnarray}&&{f}_{\mathrm{rms}}=\displaystyle \frac{{({{\rm{\Sigma }}}_{j=1}^{N}{({r}_{j}-\bar{r})}^{2}/N)}^{1/2}}{\bar{r}},\end{eqnarray} \tag{ 1 }$$

where r_j is the count rate in phase bin j and $\bar{r}$ is the phase-averaged count rate. The time history of the pulsed fraction for the six energy bands is shown in Figure 4 (left panel). The time evolution of f_rms in the individual energy bands is similar to the full-band f_rms. However, the peak pulsed fraction increases substantially with energy. The relationship between the full 0.2–12 keV band f_rms and the 0.1–10 keV luminosity (see Table 1) is shown in Figure 4 (right panel). A clear turnaround is present in both the rising and fading portions of the outburst.

Zoom In Zoom Out Reset image size

We analyzed variability on timescales other than the neutron star spin period by estimating the power spectral density (PSD; P(f)) for each epoch. First, we binned the data in 100 ms bins using photon times that had been corrected to the solar system barycenter. We computed deadtime and exposure, _i, as above, allowing an estimate of the true count rate, r_i = c_i/_i.

Because NICER only views the source for a fraction of each ISS orbit, the window function for a given epoch is typically complicated, leading to a convolution of the true power spectrum with a series of sinc functions. In the PSD, this is squared, so a pure tone will be spread into side lobes whose magnitude decreases as 1/f², a phenomenon known as spectral leakage. It is therefore critical both to eliminate as much pulsed power from the signal as possible and to pay careful attention to the resulting spectra for evidence of such leakage.

To reduce the pulsed signal, we first estimate the spin period and period derivative for each epoch, taking into account binary motion, and use this to assign a pulse phase to each time bin, ϕ(t_i). We then approximate the power at the spin frequency with a truncated Fourier series whose coefficients we estimate as, for example, ${\sum }_{i}{r}_{i}\cos (2\pi k{\phi }_{i})/{\sum }_{i}{r}_{i}$, with k between 1 and 10. We evaluate the Fourier series at each time bin using ϕ(t_i) and _i and subtract this from the observed rate.

To compute a PSD for visualization for each epoch, we selected contiguous segments of data with at least 2048 samples and then zero-padded the (mean subtracted) data to 2¹⁴ samples. Most segments were shorter than this. We computed the discrete Fourier transform using the FFT algorithm and co-added the resulting spectra using the interval lengths as weights. We normalized the spectra according to the convention common in the literature (e.g., Nowak et al. 1999), in which the PSD integrated over frequency yields the rate-normalized variance, ${f}_{\mathrm{rms}}^{2}$ (see Equation (1)). Using shorter, contiguous segments degrades spectral resolution and precludes examination of frequencies <1 mHz, but it provides superior window functions and is more robust against any uncorrected deadtime effects that evolve over an epoch.

The resulting PSDs for 16 epochs appear in Figure 5, which we explain in detail here. For display purposes, all spectra have been logarithmically smoothed with a resolution of about 0.05 dB. The tan trace in each panel of Figure 5, with an obvious harmonic series, is the PSD of the mean pulsed signal, computed as described above, multiplied by the exposure for the particular observation. The fundamental appears at the spin frequency 0.1014 Hz. It is typically quite strong compared to the broadband variations. Because the exposure (window function) is the same as for the real data, this PSD is useful for assessing the effect of the window function on pure tones. We also include a dashed gray line whose slope is fixed to 1/f² to aid in assessing spectral indices and leakage. The blue trace shows the PSD of the observed count rate after subtracting the mean pulsed component. Some remaining pulsed power is evident (e.g., in the observation on MJD 58031.0), but such power is small compared to the broadband features. The dark horizontal line shows the white noise level, which for Poisson statistics is given by $\langle r/\epsilon \rangle$. We find this to be in good agreement with the observed spectra. Finally, the black trace shows the power spectrum with white noise subtracted. This is the best empirical estimate of the aperiodic PSD and can be used for visual characterization. We note here that we also computed PSDs by omitting each MPU slice in turn, and found no substantial difference. Finally, we note that the relative power levels are high, such that any contributions from, for example, pointing jitter, are negligible by comparison.

Zoom In Zoom Out Reset image size

We have also fit functional forms to the power spectra using a maximum likelihood approach. For this purpose, we did not co-add the segments within an epoch, but rather computed FFTs with lengths exactly equal to the segments, such that the resulting FFT bins are statistically independent. The elements of the PSD then follow a ${\chi }_{2}^{2}$ distribution when scaled appropriately, so the log likelihood becomes $\mathrm{log}{ \mathcal L }=-{\sum }_{i}\mathrm{log}{P}_{m}({f}_{i},\lambda )+P({f}_{i})/{P}_{m}({f}_{i},\lambda )$, where P_m(f) is the model PSD for parameters λ, and we maximize this likelihood evaluated over the unsmoothed power spectra. We exclude spectral bins with f > 2 Hz to minimize the effects of aliased power from f > 5 Hz (the Nyquist frequency). In the literature, lorentzians of various widths are widely used to model PSDs (e.g., Reig & Nespoli 2013), but they asymptotically follow 1/f². The observed PSDs are typically shallower than this, so we have instead fit the broadband power with a power law with a low frequency cutoff, $P{(f)\propto {[1+(f/{f}_{c})}^{2}]}^{\alpha /2}$. We model narrower features with a lognormal function, and we include the white noise contribution given above. The resulting power law fits are shown with red traces, individual components appearing as dashed lines and the total model solid.

There is clear evolution in the power spectra. The earliest observation is essentially white at low frequencies with a low quality quasi-periodic oscillation (QPO)-like feature at 50 mHz. QPOs have been observed already in a number of accreting pulsars ranging from ∼1 mHz to ∼40 Hz (see, e.g., Paul & Naik 2011; Reig & Nespoli 2013), and are generally associated with inhomogeneities in the inner accretion disk. The spectral break is almost exactly at the neutron star spin frequency (0.1014 Hz), and the high-frequency behavior is slightly shallower than 1/f². Inspection of the window function shows this steeper tail cannot simply be spectral leakage from the QPO, so it reflects a distinct physical process. On MJD 58030.9, the QPO is stronger, narrower, and at higher frequency, about 70 mHz. However, because of the low spectral resolution at these frequencies, we are unable to quantitatively characterize any evolution of the QPO frequency and shape with luminosity.

As the luminosity increases, the power law cutoff disappears and the power law becomes harder with typical indices of −1.1 to −1.3. At higher frequencies, a broad peak of power, about a decade wide and largely overlapping the band including the spin frequency and its harmonics, appears (MJD 58031.0–58049.0). At these epochs, this broad peak can be identified with the lognormal component (dashed red lines) of the model fit. Interestingly, this component is not necessarily new, as the normalized power at ∼0.5 Hz is almost constant at all epochs; thus this component may instead model a deficit in power around 0.1 Hz. In the brightest part of the outburst (MJD 58049.0–58073.1), the spectra are consistent with a single power law. As the luminosity decays, the broad high-frequency peak re-appears (MJD 58089.4–58151.3). Finally, the spectrum returns to a flat low frequency pedestal terminated by a cutoff/QPO near the spin frequency and a steep (∼−1.7 to −2.0) high-frequency tail.

In many accreting X-ray pulsars, the source spectrum is found to be flux-dependent. When the observed spectrum is fit using a power-law component, the spectral photon index Γ proportionally or inversely correlates with the source luminosity, according to the accretion regime. A critical luminosity L_c distinguishes between two accretion regimes (Basko & Sunyaev 1976): at relatively low luminosities L < L_c, the height of the accretion structure is governed by the ram pressure of the infalling material, called "sub-critical," while at higher luminosities L > L_c a radiation-dominated shock rises in the magnetic column, governed by radiation pressure in the accreting material, called "super-critical." For sub-critical sources, the spectrum of many accreting X-ray pulsars has been found to harden as the observed X-ray flux increases (Klochkov et al. 2011; Fürst et al. 2014; Malacaria et al. 2015; Postnov et al. 2015; Epili et al. 2017). On the other hand, super-critical sources show the opposite behavior—that is, a spectral softening as the flux increases (Reig & Nespoli 2013; Postnov et al. 2015; Epili et al. 2017). Therefore, thanks to the opposite observed behaviors, the flux dependence of the spectral photon index can be used to identify the accretion regime.

A model-independent way to probe whether a source is accreting in the sub- or super-critical regime is through the analysis of the hardness ratio. Reig and Nespoli (2013) studied the patterns of variability of the hardness ratio during outbursts in a number of Be/X-ray pulsars. Those authors identified two branches, the horizontal branch (HB) and the diagonal branch (DB), in the hardness-intensity diagrams for these sources. The spectral hardness follows a certain path in the hardness-intensity diagram and undergoes state transitions from one branch to another—that is, it "turns over" if the source luminosity passes through a certain critical value. These diagrams used a intensity in the 4–30 keV band, using data from the RXTE Proportional Counter Array (PCA, Jahoda et al. 1996) and a soft color defined as the ratio between the count rates in two energy bands—namely, 7–10 keV and 4–7 keV.

We investigated the spectral behavior of Swift J0243.6+6124 by examining the hardness ratio evolution throughout the entire outburst, with both NICER and Fermi GBM data. NICER data are useful for reproducing the original Reig and Nespoli (2013) soft color (SC, 7–10 keV/4–7 keV), while GBM enables one to explore the hardness ratio behavior at higher energy and offers uniformly sampled coverage of the entire outburst. The 4–10 keV NICER count rate has been used as a proxy for the 4–30 keV total intensity originally employed by Reig and Nespoli (2013). For GBM, the hardness ratio has been chosen as the ratio between the pulsed flux levels in the 12–16 and 8–12 keV energy ranges, whose edges are defined by the channel edges of CSPEC data (see Section 2.2), while the chosen bands are such that the intervals are insensitive to interstellar absorption while guaranteeing sufficient statistics. Likewise, the 8–16 keV pulsed flux has been used to represent the GBM outburst intensity. Figure 6 shows the evolution of the hardness ratios with intensity for both the NICER and the GBM energy bands. GBM data clearly show both the HB at lower intensity and the DB at higher intensity, with a turnover taking place in between. However, NICER seems to trace the DB more clearly, while only a hint of the HB appears as the luminosity decreases along the decay stage of the outburst. Moreover, the turnover traced in the hardness-intensity diagram by NICER data seems to happen at a lower luminosity than GBM's. We will discuss possible implications of these results in Section 5.2.

Zoom In Zoom Out Reset image size

The energy spectra of accretion-powered X-ray pulsars can be described by simple phenomenological models despite the complex emission processes acting close to the neutron star. In support of our timing studies and to explore the nature of X-ray emission, we attempted to understand the spectral behavior of Swift J0243.6+6124 during the outburst. Photon energy spectra were extracted from all available data (see Table 1). We followed the standard analysis method as described in Sections 2.1 and 3.1. In addition to conventional filtering criteria, a SUNSHINE flag was also added to GTI that can separate the data at times when NICER was exposed to direct sunlight ("day" orbits), and times when NICER was shadowed by the Earth ("night" orbits). This was done as a first step to address the known instrument artifacts from O, Si, and Au in the NICER spectra at energies close to 0.5, 1.8, and 2.2 keV, respectively (see also Ludlam et al. 2018). These residuals were more prominent in high count rate observations, complicating spectral analysis at present for brighter sources.

We adopted a normalization technique that uses the Crab Nebula's observed spectrum to mitigate these instrument features. For this, we extracted Crab spectra from (ObsIDs 1011010101, 1011010201, and 1013010101−1013010122), in a manner similar to Swift J0243.6+6124. The products from day and night orbits were combined into respective groups using the FTOOL addspec. We accumulated time-averaged spectra of the Crab with an effective exposure of ∼26 ks for day and ∼14 ks for night orbits. Together with a background data set (RXTE blank sky field 6; Jahoda et al. 2006), and the instrument response (energy redistribution matrix and effective area files, version 0.06), the Crab spectrum was described in the 0.25–10 keV range by an absorbed power law. A poor fit was obtained at this point due to residuals in the soft X-rays arising from the preliminary instrument calibration at this stage. We froze the absorption column density from the above fit and refitted the continuum in a flat, featureless (4–10 keV) range. At these energies, the model yielded a better fit with a photon index close to 2. We next extrapolated the energy spectrum back down to 0.25 keV and created a fake spectrum using the fakeit command in XSPEC for the Crab exposure. By construction, simulated data was free of significant instrumental bias below 4 keV, and taking its ratio with the Crab spectrum, therefore, generated a "template" that carried the artifacts. Finally, we divided the Swift J0243.6+6124 spectrum by the template, using the FTOOL MATHPHA.

The 3–79 keV energy spectrum of the pulsar with NuSTAR can be expressed by a cutoff power-law, high energy cutoff power-law, and/or negative and positive exponential cutoff power-law models modified with a blackbody component (Jaisawal et al. 2018). For the purpose of estimating the pulsar's unabsorbed flux from NICER spectral analysis, we instead considered a simple cutoff power law to fit the 0.3–10 keV continuum. We noted that spectra at brighter phases of the outburst (>3 Crab intensity) suffered from instrument features, despite the Crab-renormalization technique, limiting the present study. Nevertheless, our results show a remarkable intensity variation and gradual spectral evolution across the outburst consistent with the accretion-regime transition scenario outlined in Section 3.4 and discussed in greater detail in Section 5. The unabsorbed flux has been approximately corrected for deadtime, as derived in Section 3.1 for each observation. The 0.1–10 keV luminosity was calculated by assuming isotropic emission and a source distance of 7 kpc (Bailer-Jones et al. 2018) and is presented in Table 1. A detailed spectral analysis of the pulsar will be presented in a future paper using refined calibration, response, and effective area files.

We have determined an orbital ephemeris (Table 2) for Swift J0243.6+6124 using pulse frequencies from GBM and the 15–50 keV count rates from Swift BAT. We binned the GBM Channel 1 (12–25 keV) CTIME data from MJD 58027–58150 to 250 ms and then fitted it to a semi-empirical background model. Next we subtracted the background model and extracted pulsed flux and frequencies for Swift J0243.6+6124 in the manner described in Finger et al. (1999) and Jenke et al. (2012). We performed a search for frequency and frequency derivative using pulse profiles folded over one day intervals. The intrinsic spin-up of an accreting HMXB pulsar, $\dot{\nu }$, is proportional to L^6/7 when accretion is mediated through a disk (Lamb et al. 1973). We applied the techniques described in (Tsygankov et al. 2016) to model this intrinsic spin-up using the Swift BAT rates as a proxy for the X-ray luminosity, because simple polynomial models failed to fully model the observed spin-up.

The BAT had considerable issues with measuring the rates for this source between MJD 58060 and 58075. This was caused by a known data overflow problem at very high count rates. When the counts exceeded the maximum value that could be stored, they wrapped and were instead stored as extremely low values.¹⁹ In order to correct the rates, only the peak rates during this period were used and the rates and errors corresponding to the times of measured frequencies were interpolated. The model emission frequencies were determined by integrating the modeled spin-up determined from the Swift BAT rates. The model frequencies were corrected for the current orbital parameters as described in Deeter et al. (1981). The model frequencies were then compared to the measured frequencies. Using the Metropolis–Hastings Algorithm (Hastings 1970), a search was performed to find the best-fit orbital parameters. The torque model based on the Swift BAT rates failed to produce a reasonable fit to the data due to the high spin-up during the peak of the outburst. Parameters that increased the index on the BAT rates at peak flux were added to the torque model: one parameter that increased the power index on the BAT rates above a threshold and, β, the threshold on the BAT rates in which the addition to the index would be applied. The new model for the spin-up rate was

$$\begin{eqnarray}&&\dot{\nu }\propto {R}^{{i}_{1}+{i}_{2}},\end{eqnarray} \tag{ 2 }$$

where R are the interpolated BAT rates, i₁ is the index when the BAT rates are below β, and i₁ + i₂ is the index when the BAT rates are above β.

The new model produced a better fit to the data, but there were still significant residuals remaining. Using the best-fit parameters to the orbit, the measured frequencies where corrected for the best-fit orbit and a detailed phase model was determined. Using the new phase model, a new search for pulsations was performed. A phase offset δϕ from the phase model was estimated for each 1 day mean pulse profile by fitting the Fourier amplitudes of the profile to a scaled and shifted template profile. The new pulsar orbital phase offsets were fit to the orbital model given by Deeter et al. (1981) using circular orbital elements an a polynomial background model to remove the remaining residuals. Only the tail of the outburst (between MJD 58098 and 58154) was used to avoid the large torque noise during the peak of the outburst. The refined orbital elements were used to further refine the phase model, and again a new search for pulsations was performed. This was repeated until there was no significant change in the best-fit orbital elements. To quote orbital element uncertainties representative of systematic effects, statistical errors in the frequency measurements were inflated to obtain a reduced chi-squared near 1.0 (46.6/43 dof). Sources of systematic error may include changes in the emission beam within the integration interval or throughout the outburst.

The frequency history, pulsed flux, and ephemeris for this source and all other sources monitored by the GBM Pulsar Monitoring team may be found at the GBM pulsar website.²⁰ See the Appendix for a detailed analysis of GBM energy dependent pulse profiles and their comparison with NICER.

The amount of data that NICER, Fermi GBM, and Swift have accumulated on the outburst of Swift J0243.6+6124 is very large, and it is not easy to understand every detail in light of our current understanding of the accretion flows onto magnetic neutron stars. That the neutron star is magnetic is clear from the fact that we see distinct pulses at a discrete frequency with clear indications of spin period changes over the course of the outburst. And further, those spin period changes appear to be correlated with the luminosity (thus the accretion rate onto the neutron star), lending confirmation to the idea that we are seeing rapid spin-up of a neutron star in a close orbit around a Be early-type star. The observation of a spin-up torque that is correlated with the luminosity is a clear indication that the accretion is mediated through an accretion disk. More quantitatively, an accretion disk will form around the pulsar if the matter in the Be star's circumstellar disk has a net angular momentum greater than the angular momentum of matter being spun around by the rotating magnetosphere. We can calculate the radius of the Be star's circumstellar disk R_disk using the relation from Hanuschik (1989) that relates it to the measured Hα equivalent width, EW = −10.3 Å, from Kouroubatzakis et al. (2017), to obtain R_disk ≈ 11R_c. Then using this value in the criterion from Klus et al. (2014; see also Shapiro & Lightman 1976; Wang 1981) along with parameters of the system (i.e., companion mass and radius of ∼16 M_⊙ and ∼7 R_⊙, respectively, from Bikmaev et al. 2017), a range of mass accretion rates $\dot{M}$ ∼ 10⁻¹¹–10⁻⁸ M_⊙ yr⁻¹, and an assumed magnetic field for the neutron star in the range B ∼ 10¹²–10¹⁴ G, we find that matter accreted onto the pulsar from the circumstellar disk of the Be star first forms an accretion disk around the pulsar.

Early in the outburst, before approximately MJD 58033, the pulse profile is complex, but dominated by a single asymmetric peak (see Figure 1 inset (a)). A remarkable change occurs near MJD 58037, in that the pulse shape transitions to a single simpler peak, both in the NICER (Figure 1 inset (b)) and GBM energy bands (See Appendix, Figure 8). Then, as the accretion luminosity takes off, a clear two-peak structure emerges (Figure 1 inset (c)) and dominates the emission pattern at all energy ranges (see also Figure 10). This two-peak flow arrangement lasts through the peak luminosity of ∼1.8 × 10³⁹ erg s⁻¹ at MJD 58064 (Figure 1 insets (d)–(f)) and as the luminosity goes down below ∼ 2 × 10³⁸ erg s⁻¹. Once MJD 58100 is past, both the NICER (Figure 1 inset (g)) and GBM (Figure 11) pulse profiles appear to go back to a single peak configuration. The pulses are very broad at this point, and it is not clear if the emission is from one pole or two poles that are essentially blended together in the light curve. Later on, past MJD 58135 the accretion configuration changes again in that now the soft NICER band stays single peaked (Figure 1 inset (h)), while the higher energy bands return to a complex structure dominated by a single asymmetric peak. At the same time, past MJD 58135 GBM gets the double-peaked pulse profile structure back (Figure 12).

Such extreme luminosity-dependent character of the pulse profiles is usually interpreted in terms of a change in the accretion regime. At lower luminosity the flow may be stopped primarily by coulomb collisions between particles as suggested by Basko and Sunyaev (1976; see also Becker et al. 2012). In this regime, the flow stopping region is located near the neutron star surface where the plasma density is very high while the emission from the stopping region may be dominated by a pencil beam component. However, once the accretion rate goes above a certain critical value, radiation pressure becomes the dominant stopping mechanism via photon scattering off of the electrons in the accretion flow. This is when the radiation-dominated shock transition begins to work its way up the magnetic column, and an accretion column emerges, while a fan-beam emission pattern escaping the accretion column walls dominates the emission.

The critical luminosity L_c that marks the transition between the coulomb dominated and radiation-pressure dominated accretion flow, that is the transition from the sub-critical to super-critical accretion regime, has been calculated by Becker et al. (2012) and Mushtukov et al. (2015a). The computation of the exact value of critical luminosity strongly depends on a number of parameters that reflect the accretion physics and are currently not fully understood or constrained, such as the geometrical parameter Λ characterizing whether the NS accretes from a wind or from a disk, the shape of the photon spectrum inside the column, and the magnetic field strength B. However, assuming canonical neutron star parameters, the standard critical luminosity is of the order of L_c ∼ 10³⁷ erg s⁻¹. This luminosity level is reached around MJD 58029 and 58166, very early and very late in the outburst, respectively (see Table 1). In fact, the above described pulse profile changes show clear transitions around a much higher luminosity, ∼10³⁸ erg s⁻¹, around MJD 58037 and 58135, suggesting a much higher than typical value for the critical luminosity. This finding is further corroborated by the spectral analysis and the study of the pulsed fraction discussed in the next sections, where additional observational evidence is gathered in support of the accretion regime transition scenario, and to further constrain the value L_c of the critical luminosity.

Figure 6 shows the hardness-intensity diagram for both the NICER and GBM energy bands throughout the entire outburst of Swift J0243.6+6124. This diagram represents a model-independent way to probe the correlation between the spectral hardness and the observed source intensity, which has already been demonstrated for a number of accreting pulsars (e.g., Reig & Nespoli 2013). However, the exact mechanism that leads to the spectral hardness dependence on flux is not yet clear. It is generally accepted that at the sub-critical regime (L < L_c), when the accretion rate grows the Coulomb layer that halts the accreting matter is pushed down into the accretion structure, where the electron gas becomes denser and hotter. Within this scenario, radiation is scattered up to higher energies via the inverse Compton effect, and a harder spectrum emerges. Super-critical sources (L > L_c) show the opposite behavior—that is, a spectral softening at higher flux levels. This can be understood in the framework of two different scenarios, both of them involving the rise of an accretion column. On the one hand, the turnover to softer spectra at higher luminosity is due to the plasma temperature decrease upwards in the column (Basko & Sunyaev 1976; Becker et al. 2012), which leads to the emission of softer photons. On the other hand, the taller column leads to a lesser fraction of radiation reflected by the neutron star atmosphere (a process that induces spectral hardening), and therefore to a softer spectrum (Postnov et al. 2015).

Even if the exact mechanisms responsible for the spectral dependence on flux are not yet fully understood, Reig and Nespoli (2013) showed that the turnover present between the HB and the DB is happening at a luminosity level comparable to the critical luminosity. Therefore, deriving the source luminosity at the turnover level is of key importance.

In the NICER energy bands, the hardness ratio shows a hint of a turnover near the beginning of the outburst, where a transition occurs from the HB to the DB (see Figure 6). The farthest-right SC (soft color; Section 3.4) for NICER occurred on MJD 58031.4. Recognizing that NICER observations are sparse and there is possible jitter in the SC value based on other sources in Reig and Nespoli (2013), we conservatively estimate that this turnover occurred sometime between MJD 58031 and MJD 58037. On the other hand, the red diamonds denoting the outburst decay in Figure 6 show that if indeed we saw the turnover during the outburst rise, we missed it on the decay. The more frequent observations obtained during the rise help constrain the decay turnover between MJD 58145 and MJD 58162. Comparing those events with our spectral results (see Section 3.5 and Table 1), we derive a luminosity for the turnover in the range ${L}_{0.1-10\mathrm{keV}}\sim (0.2\mbox{--}1.1)\,\times {10}^{38}$ erg s⁻¹. We note that Jaisawal et al. (2018) measured ${L}_{3-70\mathrm{keV}}\sim 5.1\times {10}^{37}$ erg s⁻¹ (corrected to d = 7 kpc) with NuSTAR on MJD 58031; thus NICER and NuSTAR derived luminosities are consistent within a factor of ∼2. As a comparison, for typical values of the neutron star parameters, Becker et al. (2012) calculate the critical luminosity as

$$\begin{eqnarray}&&{L}_{{\rm{c}}}\sim 1.5\times {10}^{37}\,{B}_{12}^{16/15}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 3 }$$

where B₁₂ is the magnetic field strength in units of 10¹² G, which is consistent with the value obtained in this work with NICER, if the magnetic field is higher than the typical ∼10¹² G.

On the other hand, Figure 6 shows that in the Fermi GBM energy bands, the turnover takes place at ∼2 keVcm⁻² s⁻¹ , corresponding to MJD 58040–58043 and 58098–58105, for the rising and the decay stages of the outburst, respectively. GBM pulsed fluxes at the turnover can be converted into a critical luminosity value by a simple scaling of the NuSTAR observed luminosity. In addition, as a consistency check, we also converted GBM pulsed fluxes into luminosity by comparing them with the Swift BAT count rates in the same time windows, and using a conversion factor of 1.54 × 10⁻⁷ erg cm⁻² cnt⁻¹ (Doroshenko et al. 2018). In the former case, the critical luminosity that results is ${L}_{3-70\mathrm{keV}}\sim (2.6\mbox{--}3.4)\times {10}^{38}\,$ erg s⁻¹ for the rising stage of the outburst, and ${L}_{3-70\mathrm{keV}}\sim (2.2\mbox{--}3.4)\times {10}^{38}\,$ erg s⁻¹ for the decay. In the latter case, the resulting critical luminosity is ${L}_{15-50\mathrm{keV}}\sim (3.2-4.3)\times {10}^{38}\,$ erg s⁻¹ for the rising stage of the outburst, and ${L}_{15-50\mathrm{keV}}\sim (3.8-4.8)\times {10}^{38}\,$ erg s⁻¹ for the decay. In any case, the turnover in the GBM energy bands happens at significantly higher luminosities than NICER's. To the best of our knowledge, we excluded any instrumental origin for such an effect, thus concluding that the dependence of the turnover luminosity on the energy band represents a physical effect. An energy- and luminosity-dependent hardness ratio could, in fact, reflect a dependence of the spectral cutoff energy on the luminosity. Previous works showed that when the spectra of accreting pulsars are fit with a cutoff power-law model, the cutoff energy is directly correlated with the power-law photon index and, at the same time, inversely/directly correlated with the luminosity for the sub-/super-critical accretion regimes, respectively (Ferrigno et al. 2013; Müller et al. 2013; Reig & Nespoli 2013; Malacaria et al. 2015). Since the cutoff energy is indicative of the electron temperature, its variability can be interpreted as the result of a Compton cooling efficiency dependence on the luminosity and the accretion regime. Furthermore, we note that the typical cutoff energy in accreting X-ray pulsars is E_cut ∼ 15–20 keV—that is, around the GBM bands used to constrain the turnover. Following this qualitative scenario, the combination of luminosity-dependent cutoff energy and photon index can lead to a difference in the turnover luminosity measured in different energy bands. However, we note that the energy bands used here to measure the turnover can be affected by other factors, such as the Fe K line complex around 6–7 keV and the so-called "10 keV feature" (Coburn et al. 2002), commonly observed among accreting X-ray pulsars.

Considering these findings, we estimate a conservative value for the critical luminosity bracketed by the values separately derived with NICER and GBM data—that is, (0.2–4.8) × 10³⁸ erg s⁻¹. Given these bounds, the turnover luminosity can be used to constrain the surface magnetic field value. Following the approximate equation for critical luminosity provided by Becker et al. (2012) and cited above, we constrain the magnetic field to lie within the range (0.2–2.6) × 10¹³ G.

In the power spectra, a QPO-like feature at 50–70 mHz is found in the 0.2–12 keV band with NICER early (MJD 58029–58031) and late (MJD 58166–58176) in the outburst. This feature can be used to estimate the magnetic field if we assume that the QPO frequency ν_QPO is equal to the Keplerian orbital frequency ν_K at distance r_K, where

$$\begin{eqnarray}&&{r}_{{\rm{K}}}={\left[\displaystyle \frac{{GM}}{{(2\pi {\nu }_{{\rm{K}}})}^{2}}\right]}^{1/3}.\end{eqnarray} \tag{ 4 }$$

Next we assume that the magnetospheric radius is equal to the Alfvén radius r_A, where

$$\begin{eqnarray}&&{r}_{{\rm{A}}}={\left(\displaystyle \frac{{\mu }^{4}}{8{GM}{\dot{M}}^{2}}\right)}^{1/7}\end{eqnarray} \tag{ 5 }$$

and μ ≈ BR³ is the magnetic dipole moment (see, e.g., Frank et al. 2002). Setting r_A = r_K yields

$$\begin{eqnarray}&&\mu =4\times {10}^{31}{\rm{G}}\,{\mathrm{cm}}^{3}{M}_{1.4}^{5/6}{\dot{M}}_{-8}^{1/2}{\left(\displaystyle \frac{{\nu }_{{\rm{K}}}}{50\mathrm{mHz}}\right)}^{-7/6},\end{eqnarray} \tag{ 6 }$$

where M_1.4 = 1.4 M_⊙ and ${\dot{M}}_{-8}={10}^{-8}\,{M}_{\odot }\,$ yr⁻¹. Finally we take R = 10 km to obtain

$$\begin{eqnarray}&&B=4\times {10}^{13}\,{\rm{G}}\ {M}_{1.4}^{-5/6}{R}_{6}^{-3}{\dot{M}}_{-8}^{1/2}{\left(\displaystyle \frac{{\nu }_{{\rm{K}}}}{50\mathrm{mHz}}\right)}^{-7/6}.\end{eqnarray} \tag{ 7 }$$

When the QPO was detected, the 0.1–10 keV luminosity was (1.2–2.0) × 10³⁷ (d/7 kpc)² for MJD 58029–58031. This corresponds to a mass accretion rate of $\dot{M}=(0.9\mbox{--}1.7)\,\times {10}^{-9}\,{M}_{\odot }\,$ yr⁻¹. Substituting this rate and the QPO frequencies of 50 and 70 mHz respectively into Equation (7) yields B ∼ 1 × 10¹³(d/7 kpc) G. Substituting this value back into Equation (3), we obtain a critical luminosity of L_c ∼ 1.7 × 10³⁸ erg s⁻¹(d/7 kpc)^16/15, consistent with and better constrained than that from the hardness ratio calculations.

Another indication of changes happening at the times corresponding to the turnover or, equivalently, to the critical luminosity is found in the analysis of the pulsed fraction evolution (Figure 4). At the earliest, rising stage of the outburst, the pulsed fraction decreases with increasing flux, down to a local minimum around MJD 58037–58042. The pulsed fraction then rises up to a maximum around the outburst peak, then decreasing again during the decay stage, down to a second local minimum at MJD 58099–58106, after which a new increase takes place, somehow symmetrically to the first half of the outburst. Above the critical luminosity, the pulse fraction clearly increases with increasing luminosity (see Figure 4, right panel); below the critical luminosity, especially as the outburst declines, the relationship is less clear. Once again, the local minimum occurrence times coincide with the turnover times, and therefore with the critical luminosity, with the pulsed fraction changes closely tracking the pulse profile changes.

The exact mechanism of the pulsed fraction dependence on the luminosity is not yet clear. To complicate the matter further, the dependence is not unique: some sources exhibit pulsed fraction decreases with increasing luminosity (e.g., 4U 0115+63, Tsygankov et al. 2007; EXO 2030+375, Epili et al. 2017), while others show the opposite trend along with a non-monotonic evolution of the pulsed fraction with luminosity (see, e.g., V 0332+53, Tsygankov et al. 2010; SMC X–2, Jaisawal & Naik 2016). Our observations suggest that the two opposing behaviors—that is, the negative and positive dependence of the pulsed fraction on luminosity—are distinctive of sub- and super-critical accretion regimes, respectively. However, we note that this is a simplistic scenario, and that more-detailed models need to be considered to take into account the diverse behavior of the pulsed fraction dependence on luminosity observed in different sources and at different accretion regimes.

ULX sources have been identified in a number of galaxies. A common definition is a luminosity >10³⁹ erg s⁻¹ (Kaaret et al. 2017). At first, some of these systems were thought to contain intermediate mass black holes, ∼100–10,000 M_⊙ (Sutton et al. 2012). In 2014, the first ULX pulsar, M82 X–2, was discovered with NuSTAR (Bachetti et al. 2014), with L_0.3–10keV = 1.8 × 10⁴⁰ erg s⁻¹ and a pulse period of 1.37 s. Since then, several other ULX pulsars have been discovered, including NGC 7793 P13, with a 0.42 s period and an observed peak luminosity of ∼10⁴⁰ erg s⁻¹ (Fürst et al. 2016; Israel et al. 2017b), and NGC 5907 ULX (Israel et al. 2017a), with a period of 1.13 s and a peak luminosity of ∼10⁴¹ erg s⁻¹. All three of these systems show rapid spin-up, have sinusoidal pulse profiles, and relatively low pulsed fractions that increase with increasing energy. They exhibit both bright and faint phases with luminosities of ∼10^40–41 erg s⁻¹ and ∼10³⁸ erg s⁻¹, respectively (Kaaret et al. 2017). All three systems were known as ULX sources before their identification as pulsars. Recently, a fourth ULX pulsar, NGC 300 ULX1, with an initial period of 32 s, rapidly spinning up, and a peak luminosity of 4.7 × 10³⁹ erg s⁻¹, was discovered (Carpano et al. 2018). This pulsar has a potential cyclotron feature at ∼13 keV, suggesting a magnetic field of B ∼ 10¹² G (Walton et al. 2018). Chandra observations show a potential proton resonance scattering feature at 4.5 keV in M51 ULX8 that implies a neutron star surface magnetic field of B ∼ 10¹⁵ G (Brightman et al. 2018), although no periodicity has been reported in this case.

With a peak 0.1–10 keV luminosity of 2 × 10³⁹ erg s⁻¹, Swift J0243.6+6124 crosses the threshold for a ULX pulsar. Incorporating the Gaia distance ranges, the dominant error on the luminosity, results in a peak luminosity of (1.2–2.6) × 10³⁹ erg s⁻¹ for the Bailer-Jones et al. (2018) ±1σ range of 5.7–8.4 kpc, and (1.5–5.6) × 10³⁹ erg s⁻¹ for the 5%–95% confidence range of 6.3–12.3 kpc (M. Ramos-Lerate 2018, private communication), all above the 10³⁹ erg s⁻¹ threshold for ULXs. Swift J0243.6+6124 also shows rapid spin-up, exhibits an increasing pulsed fraction with energy, and has a relatively short spin period of ∼9.8 s. Mushtukov et al. (2017) describe an optically thick envelope around accreting ULX pulsars that produces a smooth pulse profile for L > 10³⁹ erg s⁻¹ with a pulsed fraction that increases with energy. This optically thick envelope masks the observer's view of the neutron star surface. The profiles from Swift J0243.6+6124 are quite smooth at high luminosities but are more complex than the simple single peaked profiles observed from the other ULX pulsars, suggesting a different viewing geometry or that Swift J0243.6+6124 may be close to the threshold for this effect, similar to SMC X-3 (Koliopanos & Vasilopoulos 2018). Mushtukov et al. (2015b) proposed that ULX sources and accreting pulsars are connected, with the brightest ULX pulsars having magnetar-like fields of ∼10¹⁴ G. They define an admissible luminosity range for accreting pulsars and ULX pulsars alike, where the maximum accretion luminosity can be approximated by ${L}_{\mathrm{acc}}\approx 0.35{B}_{12}^{3/4}\times {10}^{39}\,$ erg s⁻¹, for 10¹³ G < B < 10¹⁵ G, shown in Figure 7 as the solid dark line. They also define a minimum luminosity (denoted by a dashed line), based on the propeller effect, which depends on the spin period of the pulsar. We demonstrate that Swift J0243.6+6124 is consistent with this picture, based on the luminosities measured with NICER (keeping in mind that the total accretion luminosity is likely at least a factor of 2 higher to incorporate flux >10 keV) and our estimate of Swift J0243.6+6124's magnetic field based on our critical luminosity estimate from the hardness ratios, pulse fraction, and pulse profiles, and the magnetic field estimate from the QPO measured in with NICER.

Zoom In Zoom Out Reset image size

Energy-dependent NICER pulse profiles were generated by selecting events in the 0.2–1, 1–2, 2–3, 3–5, 5–8, and 8–12 keV bands. Figures 8–12 (top panels) show the evolution of the pulse profiles with time and with energy. Error bars are standard 1-σ errors assuming Poisson statistics and are smaller than the line thickness in the brighter observations. Profiles are not phase connected. At lower average count rates (<1000 cps), early (Figure 8) and late (Figure 12) in the outburst, considerable energy-dependent pulse shape variations are observed and the pulse profile is very complex, with the largest peak closest to the deepest minimum. Then a transition occurs at 500–1000 cps (Figure 8 4th–5th column, Figure 12 2nd–5th column), where the pulse profile becomes primarily a single asymmetric peak and seems to reverse in phase as to which peak is dominant. Above about 3000 cps (last two columns of Figures 8 and 10), the profile becomes more symmetric and gradually splits into two equal peaks, and the energy dependence becomes less dramatic. The minimum deepens between the two peaks as the intensity increases. As the outburst approaches its peak >16,000 cps (Figure 9), the profile again becomes asymmetric, with the peak closer to the deeper minimum again dominating. As the outburst decays, the profiles evolve very similarly through the shapes and complexities observed during the outburst rise. This analysis demonstrates the wealth of observations made possible by NICER for an extremely bright source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Pulse profiles were generated in nine energy bands from 5 to 100 keV using GBM CTTE data. For each day of data, harmonic expansions including 24 sine and cosine terms in a pulse phase model were fitted to the GBM count rates. Times were barycentered and corrected for the orbital ephemeris described in the previous section. These harmonic expansions were converted to a more typical pulse profile and are plotted in Figures 8–12, bottom panels. Error bars are plotted for every ∼4 bins, corresponding to the independent harmonics. Each profile has been shifted so that the minimum in the 8–12 keV band is at phase zero. Comparing these profiles to the NICER profiles in the same columns, the higher energy profiles also evolve from complex at first to asymmetric and single peaked but narrower, to double peaked but more asymmetric than at lower energies, back to single peaked, and complex. In the 0.2–1.0 keV band, two distinct peaks emerge very quickly, but in the higher energy range of 8–12 keV, the emission is still dominated by one primary peak. At even higher energy, in the GBM domain the pulse profile shows the distinct peaks already at early stages (see Figure 8). When the profiles are the most complex, the phasing of the GBM and NICER profiles may not be consistent.