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

Colleen A. Wilson-Hodge, Christian Malacaria, Peter A. Jenke, Gaurava K.

Jaisawal, Matthew Kerr, Michael T. Wolff, Zaven Arzoumanian, Deepto

Chakrabarty, John P. Doty, Keith C. Gendreau, et al.

► To cite this version:

Colleen A. Wilson-Hodge, Christian Malacaria, Peter A. Jenke, Gaurava K. Jaisawal, Matthew Kerr, et al.. NICER and Fermi GBM Observations of the First Galactic Ultraluminous X-ray Pulsar Swift J0243.6+6124. The Astrophysical Journal, 2018, 863 (1), pp.9. 10.3847/1538-4357/aace60. hal-01872296

HAL Id: hal-01872296

https://hal.science/hal-01872296v1

Submitted on 19 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Colleen A. Wilson-Hodge<fnref n="1" />, Christian Malacaria<fnref n="1" />,<fnref n="2" />,<fnref n="16" />, Peter A. Jenke<fnref n="3" />, Gaurava K. Jaisawal<fnref n="4" />, Matthew Kerr<fnref n="5" />, Michael T. Wolff<fnref n="5" />, Zaven Arzoumanian<fnref n="6" />, Deepto Chakrabarty<fnref n="7" />, John P. Doty<fnref n="8" />, Keith C. Gendreau<fnref n="6" />, Sebastien Guillot<fnref n="9" />,<fnref n="10" />, Wynn C. G. Ho<fnref n="11" />,<fnref n="12" />, Beverly LaMarr<fnref n="7" />, Craig B. Markwardt<fnref n="6" />, Feryal Özel<fnref n="13" />, Gregory Y. Prigozhin<fnref n="7" />, Paul S. Ray<fnref n="5" />, Mercedes Ramos-Lerate<fnref n="14" />, Ronald A. Remillard<fnref n="7" />, Tod E. Strohmayer<fnref n="6" />, Michael L. Vezie<fnref n="7" />, and Kent S. Wood<fnref n="15" />
on behalf of the NICER Science Team

*Received 2018 April 23; revised 2018 June 18; accepted 2018 June 19; published 2018 August 6*

## Abstract

Swift J0243.6+6124 is a newly discovered Galactic Be/X-ray binary, revealed in late 2017 September in a giant outburst with a peak luminosity of $2 \times 10^{39} (d/7 \text{ kpc})^2 \text{ erg s}^{-1}$ (0.1–10 keV), with no formerly reported activity. At this luminosity, Swift J0243.6+6124 is the first known galactic ultraluminous X-ray pulsar. We describe *Neutron star Interior Composition Explorer (NICER)* and *Fermi Gamma-ray Burst Monitor (GBM)* timing and spectral analyses for this source. A new orbital ephemeris is obtained for the binary system using spin frequencies measured with GBM and 15–50 keV fluxes measured with the *Neil Gehrels Swift Observatory* Burst Alert Telescope to model the system’s intrinsic spin-up. Power spectra measured with NICER show considerable evolution with luminosity, including a quasi-periodic oscillation near 50 mHz that is omnipresent at low luminosity and has an evolving central frequency. Pulse profiles measured over the combined 0.2–100 keV range show complex evolution that is both luminosity and energy dependent. Near the critical luminosity of $L \sim 10^{38} \text{ erg s}^{-1}$, the pulse profiles transition from single peaked to double peaked, the pulsed fraction reaches a minimum in all energy bands, and the hardness ratios in both NICER and GBM show a turnover to softening as the intensity increases. This behavior repeats as the outburst rises and fades, indicating two distinct accretion regimes. These two regimes are suggestive of the accretion structure on the neutron star surface transitioning from a Coulomb collisional stopping mechanism at lower luminosities to a radiation-dominated stopping mechanism at higher luminosities. This is the highest observed (to date) value of the critical luminosity, suggesting a magnetic field of $B \sim 10^{13} \text{ G}$.

**Key words:** accretion, accretion disks – pulsars: individual (SWIFT J0243.6+6124) – X-rays: binaries

## 1. Introduction

Accretion-powered pulsars are X-ray binaries consisting of a highly magnetized neutron star (~10¹<fnref n="2" /> G) accreting from an ordinary donor star (Bildsten et al. 1997). The long-term evolution of these systems is largely set by the mass of the donor (see, e.g., Tauris & van den Heuvel 2006). Among high-mass X-ray binaries (HMXBs; donor mass $\gtrsim 8 M_\odot$), we may identify two distinct types of systems: persistent or quasi-persistent X-ray sources accreting from an OB-supergiant donor in a circular or low-eccentricity binary, and episodic X-ray transients accreting from an Oe/Be-type main sequence donor in an eccentric binary. In the latter group are the so-called Be/X-ray binaries. Two types of outbursts (Stella et al. 1986) are generally observed from these systems (e.g., Paul &

Naik 2011; Laplace et al. 2017): Type I or normal outbursts, which generally begin around the time of the neutron star’s periastron passage, and Type II or giant outbursts, which are rarer, longer, reach higher luminosities $L \ge 10^{37} \text{ erg s}^{-1}$, and can start at any orbital phase. Until recently the giant outbursts were often considered to be the result of an enlarged decretion disk around the Be star, but recent work (Monageng et al. 2017) shows no correlation between the Be disk radii and the occurrence of giant outbursts. However, the normal outbursts were found to occur generally when the Be disk truncation radius was close to or larger than the Roche critical lobe radius at periastron passage. In this paper, we describe observations of a new Be/X-ray transient pulsar that is particularly luminous, making it a superb target for detailed study.

The X-ray source Swift J0243.6+6124 was initially identified as a new transient by the *Swift Burst Alert Telescope*

<footnote n="1">¹ST12 Astrophysics Branch, NASA Marshall Space Flight Center, Huntsville, AL 35812, USA</footnote>

<footnote n="2">²Universities Space Research Association, NSSTC, Huntsville, AL 35805, USA</footnote>

<footnote n="3">³University of Alabama in Huntsville, Huntsville, AL 35805, USA</footnote>

<footnote n="4">⁴National Space Institute, Technical University of Denmark, Elektrovej, DK-2800 Lyngby, Denmark</footnote>

<footnote n="5">⁵Space Science Division, U.S. Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="6">⁶Code 662, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="7">⁷MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="8">⁸Noosi Aerospace Ltd., Billerica, MA 01821, USA</footnote>

<footnote n="9">⁹CNRS, IRAP, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="10">¹⁰Université de Toulouse, CNES, UPS-OMP, F-31028 Toulouse, France</footnote>

<footnote n="11">¹¹Department of Physics and Astronomy, Haverford College, Haverford, PA 19041, USA</footnote>

<footnote n="12">¹²Mathematical Sciences, Physics and Astronomy, and STAG Research Centre, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="13">¹³Department of Astronomy, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="14">¹⁴Vitrociset Belgium for ESA/ESAC, Camino Bajo del Castillo s/n, E-28691 Villanueva de la Caada, Spain</footnote>

<footnote n="15">¹⁵Praxis Inc., Arlington, VA 22202, USA</footnote>

<footnote n="16">¹⁶NASA Postdoctoral Fellow.</footnote>
(BAT; 15–50 keV) on 2017 October 3 (Cenko et al. 2017), with 0.2–10 keV pulsations with a 9.86 s period also detected with Swift X-ray Telescope (XRT; Kennea et al. 2017). The onset of the X-ray outburst was actually detected in the 2–10 keV band a few days earlier by the Monitor of All-sky X-ray Image (MAXI) sky monitor on September 29, but was initially misidentified as the nearby source LSI +61° 303 (Sugita et al. 2017b, 2017a), only 24' away. The X-ray pulsations were confirmed in data from the Fermi Gamma-ray Burst Monitor (GBM, 12–25 keV; Jenke & Wilson-Hodge 2017), Swift XRT (Beardmore et al. 2017), and NuSTAR (3–79 keV; Bahramian et al. 2017; Jaisawal et al. 2018). NuSTAR observations also showed that the pulse shape is strongly energy dependent in the 3–79 keV band, but no cyclotron absorption features have yet been detected in the hard X-ray spectrum (Jaisawal et al. 2018). Preliminary solutions for a 26–28 day eccentric ($e \approx 0.1$) pulsar orbit were reported by Ge et al. (2017) and Doroshenko et al. (2018), with an improved solution given by Jenke et al. (2018; see also Section 4.1 below).

A variable 6 GHz radio counterpart was also detected, with a flux of <27 µJy on October 10 and 76 ± 7 µJy on November 8 (van den Eijnden et al. 2017a, 2017b). An optical counterpart with magnitude $B = 13$, USNO-B1.0 1514+0083050, was identified by positional coincidence with the Swift XRT source (Kennea et al. 2017). Optical spectroscopy corresponds to a late Oe-type or early Be-type star (Bikmaev et al. 2017; Kouroubatzakis et al. 2017). The optical counterpart shows evidence for long-term (~1000 day) V-band variability of order 0.15 mag (Kochanek et al. 2017; Nesci 2017; Stanek et al. 2017). The distance to the source was initially estimated to be 2.5 ± 0.5 kpc, based on the estimated absolute magnitude of the stellar companion (Bikmaev et al. 2017).

However, on 2018 April 25, Gaia Data Release 2 was made available, including an online catalog<fnref n="17" /> with derived distances for a large number of sources, including USNO-B1.0 1514+0083050, with a measured parallax of 0.095 ± 0.030 mas. This catalog was based on Bailer-Jones et al. (2018), where the Gaia team applied best practices to derive distances for all of their stars. USNO-B1.0 1514+0083050 corresponds to source id 465628193526364416. The goal of the Bailer-Jones et al. (2018) study is to provide purely geometric distance estimates, independent of assumptions about the physical properties of or the interstellar extinction toward stars. The distance for Swift J0243.6+6124 from this catalog is 6.8 kpc, with a 1-σ range of 5.7–8.4 kpc. A separate analysis, performed by a member of the Gaia team, using different priors, resulted in a distance of 8 kpc, with a 5%–95% confidence range of (6.3–12.3 kpc; M. Ramos-Lerate 2018, private communication). Since these distances are similar and the ~7 kpc distance is publicly available, for purposes of this paper, we adopt a distance of 7 kpc based on the Gaia measurements.

Transient systems similar to Swift J0243.6+6124 are a particularly valuable laboratory for the study of magnetically channeled accretion, as they trace through a large dynamic range of luminosity and mass accretion rate over the course of an outburst on timescales of days to weeks, before subsiding into X-ray quiescence. The goal of this paper is to use detailed analyses of the pulse profiles, pulsed fraction, power spectra, and hardness ratios over a broad energy range to understand the

changes in accretion column geometry and characteristics as the mass accretion rate evolves in the Swift J0243.6+6124 system. In this paper, we describe 0.2–12 keV Neutron star Interior Composition Explorer (NICER) and 8–100 keV GBM observations of Swift J0243.6+6124 obtained throughout the full outburst. Swift J0243.6+6124 evolved over nearly three orders of magnitude in luminosity during this outburst, reaching a peak count rate in the NICER energy band of 63,000 cts s⁻¹ (see Figure 1), nearly six times the nominal flux from the Crab Nebula and pulsar (11,000 cts s⁻¹). As the intensity increased, the pulsar underwent dramatic spin-up, as measured with Fermi GBM (see Figure 2). We describe the extremely detailed 0.2–12 keV (Section 3.2) and 8–100 keV pulse profile data sets obtained by NICER and GBM, respectively, which show considerable evolution with luminosity and energy. Detailed plots of pulse profiles versus energy as a function of time are shown in the Appendix. Like the pulse profiles, the pulsed fraction (Section 3.2) and power spectra (Section 3.3) measured with NICER also evolve with luminosity. Hardness ratio histories measured with NICER and GBM (Section 3.4) show a turnover, above which hardness is anticorrelated with luminosity. Luminosities in the 0.1–10 keV band were estimated from preliminary spectral fits to the NICER data (Section 3.5). Implications of these results, which suggest a transition in accretion regimes, are discussed in Section 5.

## 2. Observations of Swift J0243.6+6124

### 2.1. NICER Instrument and Data

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. TheXTI 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 theXTI 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.

### 2.2. Fermi GBM Instrument and Data

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

<footnote n="17">¹⁷ http://gaia.ari.uni-heidelberg.de/tap.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Outburst light curve and pulse profile evolution for Swift J0243.6+6124. The large plot shows the evolution of the average count rate for each NICER observation in the 0.2–12 keV range. The insets to the left and right show the dramatic 0.2–12 keV pulse shape variations as the outburst progresses. For the inset plots, time goes from the bottom to the top of the page for the rising portion of the outburst (left) and down from the top to the bottom of the page for the declining part of the outburst (right).</figcaption></figure>

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.

## 2.3. Swift BAT Instrument and Data

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<fnref n="18" /> (15–50 keV), provided by the BAT team, to model the accretion torques applied to the neutron star.

# 3. NICER Observation Analysis

## 3.1. NICER Event Deadtime Analysis

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 <0°015 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,

<footnote n="18">¹⁸ http://heasarc.gsfc.nasa.gov/docs/swift/results/transients/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: barycentered and orbit-corrected spin-frequency history measured with GBM. Bottom: pulsed flux measured with GBM in nine energy bands. Swift J0243.6+6124 was detected from 2017 October 1 through 2018 February 25 (MJD 58027–58174) in single day integrations and through 2018 March 3 (MJD 58180) in 3 day integrations.</figcaption></figure>

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.

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) \times 16/(25.8 \times 10^6 \text{ 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.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Swift J0243.6+6124 folded pulse profile (two cycles plotted) of one of the brightest NICER observations, ObsID # 1050390115, on 2017 November 7. The black histogram and left-hand y-axis denote the folded pulse profile in the 0.2–12 keV band, averaged over the 7 MPU slices, after corrections for GTI exposure and deadtime have been applied. The green curve and right-hand y-axis denote the deadtime corrected exposure per phase bin, averaged over the 7 MPU slices for this observation. The total deadtime corrected exposure, averaged over MPU slices, for this observation is 1276.5 s. If the bin-to-bin variations were not present, each bin would have an exposure of 12.765 s. The exposure per bin differs for each observation and shows bin-to-bin variations >1% compared to uniform exposure for total count rates >2000 counts s⁻¹.</figcaption></figure>

Fourth, to obtain the deadtime corrected exposure time for each phase bin, $\epsilon$, 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.

## 3.2. NICER Pulse Profile Analysis

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

$$f_{rms} = \frac{(\sum_{j=1}^{N}(r_j - \bar{r})^2/N)^{1/2}}{\bar{r}}, \quad (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.

## 3.3. NICER Power Spectra

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, $\epsilon_i$, as above, allowing an estimate of the true count rate, $r_i = c_i/\epsilon_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^2$, 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, $\phi(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
<table><caption>**Table 1**
NICER Observations of Swift J0243.6+6124</caption><thead><tr><th>NICER ObsID</th><th>Date</th><th>MJD Epoch</th><th>On-source Time (s)</th><th>Deadtime Corrected Average Exposure Per MPU Slice (s)</th><th>Average Rate 0.2–12 keV (counts s<sup>-1</sup>)</th><th>Luminosity<sup>a</sup> 0.1–10 keV (erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1050390101</td><td>2017 Oct 03</td><td>58029.805</td><td>3099</td><td>2643.4</td><td>164.8</td><td>1.19E+37</td></tr><tr><td>1050390102</td><td>2017 Oct 04</td><td>58030.922</td><td>4171</td><td>2782.6</td><td>239.9</td><td>1.75E+37</td></tr><tr><td>1050390103</td><td>2017 Oct 04</td><td>58031.434</td><td>21590</td><td>16645.4</td><td>281.9</td><td>2.05E+37</td></tr><tr><td>1050390104</td><td>2017 Oct 06</td><td>58032.456</td><td>24029</td><td>19320.5</td><td>423.7</td><td>2.99E+37</td></tr><tr><td>1050390105</td><td>2017 Oct 11</td><td>58037.128</td><td>2497</td><td>2198.3</td><td>1782.9</td><td>1.08E+38</td></tr><tr><td>1050390106</td><td>2017 Oct 12</td><td>...</td><td>41</td><td>0</td><td>...</td><td>...</td></tr><tr><td>1050390107</td><td>2017 Oct 16</td><td>58042.455</td><td>1495</td><td>420.9</td><td>3686.8</td><td>2.01E+38</td></tr><tr><td>1050390108</td><td>2017 Oct 21</td><td>58047.663</td><td>1394</td><td>615.7</td><td>4979.0</td><td>2.86E+38</td></tr><tr><td>1050390109</td><td>2017 Oct 22</td><td>...</td><td>2334</td><td>0</td><td>...</td><td>...</td></tr><tr><td>1050390110</td><td>2017 Oct 23</td><td>58049.271</td><td>900</td><td>735.4</td><td>6312.7</td><td>3.31E+38</td></tr><tr><td>1050390111</td><td>2017 Oct 24</td><td>...</td><td>689</td><td>0</td><td>...</td><td>...</td></tr><tr><td>1050390112</td><td>2017 Oct 25</td><td>58051.364</td><td>1003</td><td>495.8</td><td>7484.5</td><td>4.01E+38</td></tr><tr><td>1050390113</td><td>2017 Nov 02</td><td>58058.386</td><td>4566</td><td>3126.0</td><td>26583.5</td><td>1.31E+39</td></tr><tr><td>1050390114</td><td>2017 Nov 06</td><td>58063.989</td><td>311</td><td>257.4</td><td>36333.8</td><td>1.77E+39</td></tr><tr><td>1050390115</td><td>2017 Nov 07</td><td>58064.414</td><td>1495</td><td>1276.5</td><td>33175.4</td><td>1.66E+39</td></tr><tr><td>1050390116</td><td>2017 Nov 08</td><td>58065.381</td><td>3816</td><td>1915.8</td><td>31833.4</td><td>1.58E+39</td></tr><tr><td>1050390117</td><td>2017 Nov 14</td><td>...</td><td>18</td><td>0</td><td>...</td><td>...</td></tr><tr><td>1050390118</td><td>2017 Nov 15</td><td>58072.846</td><td>405</td><td>374.9</td><td>20774.0</td><td>1.32E+39</td></tr><tr><td>1050390119</td><td>2017 Nov 16</td><td>58073.460</td><td>3917</td><td>2171.9</td><td>18636.9</td><td>9.98E+38</td></tr><tr><td>1050390120</td><td>2017 Nov 17</td><td>58074.263</td><td>3916</td><td>2240.9</td><td>17297.2</td><td>1.03E+39</td></tr><tr><td>1050390121</td><td>2017 Nov 17</td><td>58074.737</td><td>1610</td><td>1517.3</td><td>16359.3</td><td>9.40E+38</td></tr><tr><td>1050390122</td><td>2017 Dec 02</td><td>58089.626</td><td>3756</td><td>3676.9</td><td>5864.4</td><td>3.25E+38</td></tr><tr><td>1050390123</td><td>2017 Dec 04</td><td>58091.203</td><td>1291</td><td>1016.0</td><td>5341.7</td><td>2.70E+38</td></tr><tr><td>1050390124</td><td>2017 Dec 05</td><td>58092.558</td><td>1704</td><td>1590.0</td><td>4943.2</td><td>2.71E+38</td></tr><tr><td>1050390125</td><td>2017 Dec 06</td><td>58093.528</td><td>1887</td><td>1854.5</td><td>4721.4</td><td>2.65E+38</td></tr><tr><td>1050390126</td><td>2017 Dec 11</td><td>58098.909</td><td>210</td><td>109.0</td><td>4187.3</td><td>2.22E+38</td></tr><tr><td>1050390127</td><td>2017 Dec 12</td><td>58099.463</td><td>733</td><td>640.9</td><td>4116.8</td><td>2.23E+38</td></tr><tr><td>1050390128</td><td>2017 Dec 13</td><td>58100.646</td><td>482</td><td>263.5</td><td>3853.4</td><td>2.21E+38</td></tr><tr><td>1050390129</td><td>2017 Dec 19</td><td>58106.402</td><td>1034</td><td>594.0</td><td>3209.5</td><td>1.84E+38</td></tr><tr><td>1050390130</td><td>2017 Dec 20</td><td>58107.913</td><td>528</td><td>279.7</td><td>3000.0</td><td>1.61E+38</td></tr><tr><td>1050390131</td><td>2017 Dec 21</td><td>58108.237</td><td>718</td><td>618.5</td><td>3110.9</td><td>1.77E+38</td></tr><tr><td>1050390132</td><td>2018 Jan 07</td><td>58125.646</td><td>210</td><td>208.1</td><td>2120.3</td><td>1.29E+38</td></tr><tr><td>1050390133</td><td>2018 Jan 08</td><td>58126.691</td><td>990</td><td>645.3</td><td>2020.6</td><td>1.30E+38</td></tr><tr><td>1050390134</td><td>2018 Jan 09</td><td>58127.622</td><td>5319</td><td>3193.7</td><td>1823.4</td><td>1.09E+38</td></tr><tr><td>1050390135</td><td>2018 Jan 10</td><td>58128.492</td><td>3133</td><td>2029.3</td><td>1844.4</td><td>1.11E+38</td></tr><tr><td>1050390136</td><td>2018 Jan 10</td><td>58129.326</td><td>3288</td><td>2686.1</td><td>1775.0</td><td>1.00E+38</td></tr><tr><td>1050390137</td><td>2018 Jan 12</td><td>58130.323</td><td>4541</td><td>3368.9</td><td>1755.2</td><td>1.04E+38</td></tr><tr><td>1050390138</td><td>2018 Jan 17</td><td>58135.850</td><td>3230</td><td>3207.7</td><td>1399.6</td><td>8.76E+37</td></tr><tr><td>1050390139</td><td>2018 Jan 27</td><td>58145.715</td><td>2266</td><td>2127.0</td><td>941.5</td><td>5.89E+37</td></tr><tr><td>1050390140</td><td>2018 Feb 02</td><td>58151.352</td><td>2114</td><td>2104.8</td><td>581.6</td><td>3.83E+37</td></tr><tr><td>1050390141</td><td>2018 Feb 13</td><td>58162.251</td><td>3636</td><td>3492.7</td><td>236.6</td><td>1.66E+37</td></tr><tr><td>1050390142</td><td>2018 Feb 17</td><td>58166.141</td><td>2768</td><td>2755.4</td><td>190.8</td><td>1.30E+37</td></tr><tr><td>1050390143</td><td>2018 Feb 20</td><td>58169.900</td><td>2326</td><td>1564.1</td><td>131.5</td><td>8.57E+36</td></tr><tr><td>1050390144</td><td>2018 Feb 26</td><td>58175.852</td><td>2175</td><td>1805.5</td><td>49.8</td><td>3.00E+36</td></tr><tr><td>1050390145</td><td>2018 Feb 27</td><td>58176.206</td><td>7058</td><td>6964.3</td><td>45.2</td><td>2.81E+36</td></tr><tr><td>1050390146</td><td>2018 Mar 04</td><td>58181.889</td><td>1535</td><td>1442.7</td><td>45.7</td><td>2.49E+36</td></tr><tr><td>1050390147</td><td>2018 Mar 25</td><td>58202.420</td><td>1251</td><td>1077.7</td><td>22.0</td><td>1.20E+36</td></tr><tr><td>1050390148</td><td>2018 Mar 27</td><td>58204.703</td><td>1705</td><td>1649.0</td><td>14.8</td><td>6.85E+35</td></tr><tr><td>1050390149</td><td>2018 Apr 09</td><td>58217.824</td><td>5003</td><td>3966.1</td><td>104.3</td><td>6.16E+36</td></tr><tr><td>1050390150</td><td>2018 Apr 10</td><td>58218.145</td><td>2025</td><td>1917.4</td><td>108.7</td><td>6.22E+36</td></tr><tr><td>1050390151</td><td>2018 Apr 11</td><td>58219.661</td><td>2278</td><td>1699.0</td><td>122.5</td><td>7.03E+36</td></tr><tr><td>1050390152</td><td>2018 Apr 12</td><td>58220.079</td><td>1606</td><td>1111.4</td><td>120.2</td><td>6.76E+36</td></tr><tr><td>1050390153</td><td>2018 Apr 13</td><td>58221.970</td><td>500</td><td>370.4</td><td>98.0</td><td>5.06E+36</td></tr><tr><td>1050390154</td><td>2018 Apr 14</td><td>58222.491</td><td>19208</td><td>16463.2</td><td>102.4</td><td>5.58E+36</td></tr><tr><td>1050390155</td><td>2018 Apr 15</td><td>58223.492</td><td>14521</td><td>12004.6</td><td>90.7</td><td>4.81E+36</td></tr><tr><td>1050390156</td><td>2018 Apr 16</td><td>58224.489</td><td>7735</td><td>6281.5</td><td>78.6</td><td>4.10E+36</td></tr><tr><td>1050390157</td><td>2018 Apr 17</td><td>58225.260</td><td>5955</td><td>5513.6</td><td>69.8</td><td>3.57E+36</td></tr><tr><td>1050390158</td><td>2018 Apr 20</td><td>58228.186</td><td>2785</td><td>2621.7</td><td>47.5</td><td>2.47E+36</td></tr><tr><td>1050390159</td><td>2018 Apr 23</td><td>58231.951</td><td>492</td><td>459.3</td><td>24.4</td><td>9.17E+35</td></tr><tr><td>1050390160</td><td>2018 Apr 24</td><td>58232.383</td><td>643</td><td>149.1</td><td>18.6</td><td>6.28E+35</td></tr></tbody></table>

**Note.**

10. We evaluate the Fourier series at each time bin using $\phi(t_i)$ and $\epsilon_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^{14}$ 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

<footnote>a Luminosity calculation assumes isotropic emission and distance = 7 kpc.</footnote>
weights. We normalized the spectra according to the conven-
tion common in the literature (e.g., Nowak et al. 1999), in
which the PSD integrated over frequency yields the rate-
normalized variance, $f_{\text{rms}}^2$ (see Equation (1)). Using shorter,
contiguous segments degrades spectral resolution and pre-
cludes 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^2$ 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.

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 $\log \mathcal{L} = -\sum_i \log P_m(f_i, \lambda) + P(f_i)/P_m(f_i, \lambda)$, where $P_m(f)$ is the model PSD for parameters $\lambda$, 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^2$. 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^2$. 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.

3.4. NICER *Hardness Ratios*

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
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (Left): Swift J0243.6+6124 root mean squared pulsed fraction in several energy bands versus time from NICER observations. (Right): Swift J0243.6+6124 0.1–10 keV luminosity ($d = 7$ kpc) versus root mean squared pulsed fraction (0.2–12 keV). The rising portion of the outburst and the fading portion are denoted by blue squares and red diamonds, respectively. During the main peak of the outburst, the pulsed fraction generally increases with increasing count rate and increasing energy, especially after a striking turnaround near $L > 10^{38}$ ($d/7$ kpc)$^2$. Error bars are smaller than the plotted symbols in both panels.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

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.

### 3.5. NICER Spectral Analysis

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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Power spectral densities for 16 epochs with NICER. The MJD and count rate (from Table 1) are given for each observation. The main features are described in the text and are, briefly, light tan, estimated pulsed PSD; blue, smoothed PSD; dark horizontal line, white noise level; black, smoothed and noise-subtracted PSD; red, total (solid) and component (dashed) fits; and dashed gray line, a $1/f^2$ line for visual comparison.</figcaption></figure>

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.

# 4. GBM Observations of Swift J0243.6+6124

## 4.1. GBM Orbital Analysis

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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Hardness ratio evolution with flux for Swift J0243.6+6124. Squares show the rise stage of the outburst, and diamonds show the decay. Both NICER (blue squares and red diamonds connected by a solid line) and Fermi GBM (green squares and cyan diamonds connected by dashed line segments) hardness ratios are shown. Error bars are smaller than the plotted symbols. Note the different scales on the y-axis: logarithmic for NICER counts, linear for GBM pulsed fluxes. Spectral evolution and hysteresis at high luminosity are evident.</figcaption></figure>

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{v}$, 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.<fnref n="19" /> 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, $\beta$, the threshold on the BAT rates in which the addition to the index would be applied. The new model for the

<table><caption>**Table 2**
Orbital Ephemeris for Swift J0243.6+6124</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>P<sub>orb</sub></td><td>27.587 ± 0.016</td><td>days</td></tr><tr><td>T<sub>π/2</sub></td><td>58115.6640 ± 0.0049</td><td>MJD</td></tr><tr><td>a<sub>s</sub>sin i</td><td>115.84 ± 0.32</td><td>lt-s</td></tr><tr><td>e</td><td>0.09848 ± 0.00042</td><td>...</td></tr><tr><td>ω</td><td>286.44 ± 0.16</td><td>degrees</td></tr><tr><td>f(M)</td><td>2.138 ± 0.015</td><td>M<sub>☉</sub></td></tr><tr><td>χ<sup>2</sup></td><td>46.6/43 dof</td><td>...</td></tr></tbody></table>

spin-up rate was

$$ \dot{v} \propto R^{i_1+i_2}, \qquad (2) $$

where *R* are the interpolated BAT rates, $i_1$ is the index when the BAT rates are below $\beta$, and $i_1 + i_2$ is the index when the BAT rates are above $\beta$.

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 $\delta\phi$ 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.<fnref n="20" /> See the Appendix for a detailed analysis of GBM energy dependent pulse profiles and their comparison with NICER.

# 5. Discussion

## 5.1. Pulse Profile Evolution at Different Accretion Regimes

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

<footnote n="19">¹⁹ https://swift.gsfc.nasa.gov/results/transients/weak/SwiftJ0243.6p6124/</footnote>

<footnote n="20">²⁰ http://gammaray.msfc.nasa.gov/gbm/science/pulsars.html</footnote>
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_{\text{disk}}$ using the relation from Hanuschik (1989) that relates it to the measured $H\alpha$ equivalent width, $EW = -10.3$ Å, from Kouroubatzakis et al. (2017), to obtain $R_{\text{disk}} \approx 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_{\odot}$ and ~7 $R_{\odot}$, respectively, from Bikmaev et al. 2017), a range of mass accretion rates $\dot{M} \sim 10^{-11}-10^{-8} M_{\odot} \text{yr}^{-1}$, and an assumed magnetic field for the neutron star in the range $B \sim 10^{12}-10^{14}$ 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^{39}$ erg s⁻¹ at MJD 58064 (Figure 1 insets (d)–(f)) and as the luminosity goes down below ~2 × $10^{38}$ 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 $\Lambda$ 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 \sim 10^{37}$ 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^{38}$ 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.

## 5.2. Hardness Ratio Evolution

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\,\text{keV}} \sim (0.2-1.1) \times 10^{38} \, \text{erg s}^{-1}$. We note that Jaisawal et al. (2018) measured $L_{3-70\,\text{keV}} \sim 5.1 \times 10^{37} \, \text{erg s}^{-1}$ (corrected to $d = 7$ kpc) with NuSTAR on MJD 58031; thus NICER and NuSTAR derived luminosities are consistent within a factor of $\sim 2$. As a comparison, for typical values of the neutron star parameters, Becker et al. (2012) calculate the critical luminosity as

$$L_c \sim 1.5 \times 10^{37} B_{12}^{16/15} \, \text{erg s}^{-1}, \quad (3)$$

where $B_{12}$ is the magnetic field strength in units of $10^{12}$ G, which is consistent with the value obtained in this work with NICER, if the magnetic field is higher than the typical $\sim 10^{12}$ G.

On the other hand, Figure 6 shows that in the Fermi GBM energy bands, the turnover takes place at $\sim 2$ keVcm$^{-2}$ s$^{-1}$, 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 \times 10^{-7} \, \text{erg cm}^{-2} \, \text{cnt}^{-1}$ (Doroshenko et al. 2018). In the former case, the critical luminosity that results is $L_{3-70\,\text{keV}} \sim (2.6-3.4) \times 10^{38} \, \text{erg s}^{-1}$ for the rising stage of the outburst, and $L_{3-70\,\text{keV}} \sim (2.2-3.4) \times 10^{38} \, \text{erg s}^{-1}$ for the decay. In the latter case, the resulting critical luminosity is $L_{15-50\,\text{keV}} \sim (3.2-4.3) \times 10^{38} \, \text{erg s}^{-1}$ for the rising stage of the outburst, and $L_{15-50\,\text{keV}} \sim (3.8-4.8) \times 10^{38} \, \text{erg s}^{-1}$ 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_{\text{cut}} \sim 15-20 \, \text{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) \times 10^{38} \, \text{erg s}^{-1}$. 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) \times 10^{13} \, \text{G}$.

### 5.3. Constraints on the Magnetic Field from the QPO-like Feature

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 $\nu_{\text{QPO}}$ is equal to the Keplerian orbital frequency $\nu_K$ at distance $r_K$, where

$$r_K = \left[ \frac{GM}{(2\pi\nu_K)^2} \right]^{1/3}. \quad (4)$$

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

$$r_A = \left( \frac{\mu^4}{8GM\dot{M}^2} \right)^{1/7} \quad (5)$$

and $\mu \approx BR^3$ is the magnetic dipole moment (see, e.g., Frank et al. 2002). Setting $r_A = r_K$ yields

$$\mu = 4 \times 10^{31} \, \text{G cm}^3 M_{1.4}^{5/6} \dot{M}_{-8}^{1/2} \left( \frac{\nu_K}{50 \, \text{mHz}} \right)^{-7/6}, \quad (6)$$

where $M_{1.4} = 1.4 M_\odot$ and $\dot{M}_{-8} = 10^{-8} M_\odot \, \text{yr}^{-1}$. Finally we take $R = 10$ km to obtain

$$B = 4 \times 10^{13} \, \text{G} \, M_{1.4}^{-5/6} R_6^{-3} \dot{M}_{-8}^{1/2} \left( \frac{\nu_K}{50 \, \text{mHz}} \right)^{-7/6}. \quad (7)$$

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

### 5.4. Pulsed Fraction Evolution

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.

## 5.5. Swift J0243.6+6124 as the First Known Galactic Ultraluminous X-Ray (ULX) Pulsar

ULX sources have been identified in a number of galaxies. A common definition is a luminosity $>10^{39}$ 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–10\text{keV}} = 1.8 \times 10^{40}$ 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^{40}$ 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^{41}$ 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^{38}$ 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 \times 10^{39}$ erg s⁻¹, was discovered (Carpano et al. 2018). This pulsar has a potential cyclotron feature at ~13 keV, suggesting a magnetic field of $B \sim 10^{12}$ 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 \sim 10^{15}$ G (Brightman et al. 2018), although no periodicity has been reported in this case.

With a peak 0.1–10 keV luminosity of $2 \times 10^{39}$ 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^{39}$ erg s⁻¹ for the Bailer-Jones et al. (2018) ±1σ range of

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Adapted from Mushtukov et al. (2015b). The solid black line represents the maximum luminosity for magnetized NSs as a function of magnetic field strength. The dotted line and shaded region above it indicates the upper limit where the spherization radius is smaller than the magnetospheric radius. Dashed lines represent the lower limit on the X-ray luminosity due to the propeller effect. The dashed-dotted line represents the critical luminosity function. The range of observed luminosity for Swift J0243 (for a distance of 7 kpc) is indicated by an orange vertical dashed line at the magnetic field value derived from the QPO, and the blue–gray shaded box indicates the estimated range of magnetic field values derived from the critical luminosity (see the text in Sections 5.2 and 5.3). For comparison, Eddington luminosity (for a 1.4 $M_⊙$ NS) is also indicated (long dash-double dotted horizontal turquoise line).</figcaption></figure>

5.7–8.4 kpc, and (1.5–5.6) × $10^{39}$ erg s⁻¹ for the 5%–95% confidence range of 6.3–12.3 kpc (M. Ramos-Lerate 2018, private communication), all above the $10^{39}$ 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^{39}$ 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^{14}$ G. They define an admissible luminosity range for accreting pulsars and ULX pulsars alike, where the maximum accretion luminosity can be approximated by $L_{\text{acc}} \approx 0.35 B_{12}^{3/4} \times 10^{39}$ erg s⁻¹, for $10^{13}$ G < $B < 10^{15}$ 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.

## 6. Summary

Swift J0243.6+6124 underwent a giant outburst lasting about 150 days from 2017 October to 2018 February, peaking well above Eddington luminosity at $1.8 \times 10^{39} \text{ erg s}^{-1}$ (0.1–10 keV). This luminosity places it in the ULX category, making it the first known galactic ULX pulsar. Strong spin-up torques, high luminosities, and inferred mass accretion rates indicate that the pulsar was accreting from a disk, not a wind.

Near a luminosity of $10^{38} \text{ erg s}^{-1}$, we identify a transition indicated by the following observations:

1. The pulse profiles (0.2–100 keV) evolve from single peaked to two distinct peaks in both NICER and GBM data.

2. The pulsed fraction reaches a minimum and then increases with increasing intensity. Increases in pulsed fraction with energy are also seen, with the fraction approaching 100% in the 8–12 keV band at the outburst peak.

3. The power spectra evolve with increasing intensity, with a QPO-like feature and a broad high-frequency peak below $L \sim 10^{38} \text{ erg s}^{-1}$, becoming consistent with a single power law in the brightest part of the outburst.

4. The hardness ratios become anticorrelated with intensity above $L \sim 10^{38} \text{ erg s}^{-1}$ in both NICER and GBM data.

All of these behaviors repeat as the source passes through the same luminosity ranges during the outburst decay. We interpret this as evidence for two accretion regimes that depend on the accretion luminosity. The first regime shows a somewhat harder spectrum in the NICER energy range and gives way to a somewhat softer spectrum as the source luminosity goes above $\sim 10^{38} \text{ erg s}^{-1}$. We have identified the transition between these two regimes as the accretion structure on the neutron star surface transitioning from a Coulomb collisional stopping mechanism to a radiation-dominated stopping mechanism, and have estimated that this occurs at a critical luminosity of $L_c = (0.2–4.8) \times 10^{38} \text{ erg s}^{-1}$, based on hardness ratios, and $L_c \sim 1.7 \times 10^{38} \text{ erg s}^{-1}$, using the QPO features to estimate the magnetic field. This lends support to models for the accretion flows put forward by Becker et al. (2012) and Mushtukov et al. (2015a). A critical luminosity of $\sim 10^{38} \text{ erg s}^{-1}$ is the highest measured to date, suggesting a higher than average magnetic field of $\sim 10^{13} \text{ G}$ for Swift J0243.6+6124.

This work was supported by NASA through the NICER mission and the Astrophysics Explorers Program. This work was also supported by NASA through the Fermi Guest Investigator Program. G.K.J. acknowledges support from the Marie Skłodowska-Curie Actions grant no. 713683 (H2020; COFUNDPostdocDTU). This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. Work at the Naval Research Laboratory by M.T.W., P.S.R., and M.K. was supported by

NASA. C.M.'s research was supported by the NASA Postdoctoral Program at the MSFC, administered by USRA.

Facilities: NICER, Fermi, Swift.

# Appendix
## Detailed Energy Dependent Pulse Profile Analysis

### A.1. NICER Pulse Profile Analysis

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-\sigma$ 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.

### A.2. GBM Pulse Profile Analysis

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.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** (Top panel): Swift J0243.6+6124 pulse profiles from NICER observations in six energy bands. Profiles are arbitrarily aligned so that the minimum of the 0.2–12 keV profile is at phase 0.0. (Bottom panel) Swift J0243.6+6124 1 day pulse profiles measured with GBM in seven energy bands. Profiles are arbitrarily shifted so that the minimum in the GBM 8–12 keV band (not shown) is at phase zero. (Both panels): the NICER and GBM profiles are color coded using average NICER 0.2–12 keV count rate for each observation. The NICER and GBM profiles are folded with the same ephemeris, using the epoch in MJD that is the title of each column. This figure includes observations from MJD 58029.8 to 58049.3, the rise of the outburst.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (Top panel): Swift J0243.6+6124 pulse profiles from NICER observations in six energy bands. Profiles are arbitrarily aligned so that the minimum of the 0.2–12 keV profile is at phase 0.0. (Bottom panel) Swift J0243.6+6124 1 day pulse profiles measured with GBM in seven energy bands. Profiles are arbitrarily shifted so that the minimum in the GBM 8–12 keV band (not shown) is at phase zero. (Both panels): the NICER and GBM profiles are color coded using average NICER 0.2–12 keV count rate for each observation. The NICER and GBM profiles are folded with the same ephemeris, using the epoch in MJD that is the title of each column. The epoch in MJD is the title of each column. This figure spans MJD 58051.3–58074.7, including the peak of the outburst.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (Top panel): Swift J0243.6+6124 pulse profiles from NICER observations in six energy bands. Profiles are arbitrarily aligned so that the minimum of the 0.2–12 keV profile is at phase 0. (Bottom panel) Swift J0243.6+6124 1 day pulse profiles measured with GBM in seven energy bands. Profiles are arbitrarily shifted so that the minimum in the GBM 8–12 keV band (not shown) is at phase zero. (Both panels): the NICER and GBM profiles are color coded using average NICER 0.2–12 keV count rate for each observation. The NICER and GBM profiles are folded with the same ephemeris, using the epoch in MJD that is the title of each column. This figure includes MJD 58074.7–58100.6, the beginning of the declining phase of the outburst.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (Top panel): Swift J0243.6+6124 pulse profiles from NICER observations in six energy bands. Profiles are arbitrarily aligned so that the minimum of the 0.2–12 keV profile is at phase 0.0. (Bottom panel) Swift J0243.6+6124 1 day pulse profiles measured with GBM in seven energy bands. Profiles are arbitrarily shifted so that the minimum in the GBM 8–12 keV band (not shown) is at phase zero. (Both panels): the NICER and GBM profiles are color coded using average NICER 0.2–12 keV count rate for each observation. The NICER and GBM profiles are folded with the same ephemeris, using the epoch in MJD that is the title of each column. This figure includes MJD 58106.4–58129.3, as the outburst continues to fade.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** (Top panel): Swift J0243.6+6124 pulse profiles from NICER observations in six energy bands. Profiles are arbitrarily aligned so that the minimum of the 0.2–12 keV profile is at phase 0.0. (Bottom panel) Swift J0243.6+6124 1 day pulse profiles measured with GBM in seven energy bands. Profiles are arbitrarily shifted so that the minimum in the GBM 8–12 keV band (not shown) is at phase zero. (Both panels): the NICER and GBM profiles are color coded using average NICER 0.2–12 keV count rate for each observation. The NICER and GBM profiles are folded with the same ephemeris, using the epoch in MJD that is the title of each column. This figure includes MJD 58130.3–58175.9 for NICER and MJD 58130.3–58166.1 for GBM, as the outburst continues to slowly fade. Swift J0243.6+6124 dropped below GBM’s one day detection threshold on MJD 58169.</figcaption></figure>
ORCID iDs

Colleen A. Wilson-Hodge https://orcid.org/0000-0002-8585-0084

Gaurava K. Jaisawal https://orcid.org/0000-0002-6789-2723

Matthew Kerr https://orcid.org/0000-0002-0893-4073

Michael T. Wolff https://orcid.org/0000-0002-4013-5650

Deepto Chakrabarty https://orcid.org/0000-0001-8804-8946

Sebastien Guillot https://orcid.org/0000-0002-6449-106X

Wynn C. G. Ho https://orcid.org/0000-0002-6089-6836

Craig B. Markwardt https://orcid.org/0000-0001-9803-3879

Paul S. Ray https://orcid.org/0000-0002-5297-5278

Ronald A. Remillard https://orcid.org/0000-0003-4815-0481

Tod E. Strohmayer https://orcid.org/0000-0001-7681-5845

References

Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Natur, 514, 202

Bahramian, A., Kennea, J. A., & Shaw, A. W. 2017, ATel, 10866

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, arXiv:1804.10121

Basko, M. M., & Sunyaev, R. A. 1976, MNRAS, 175, 395

Beardmore, A. P., Breeveld, A. A., Gropp, J. D., et al. 2017, GCN, 22020

Becker, P. A., Klochkov, D., Schönherr, G., et al. 2012, A&A, 544, A123

Bikmaev, I., Shimansky, V., Irtuganov, E., et al. 2017, ATel, 10968

Bildsten, L., Chakrabarty, D., Chiu, J., et al. 1997, ApJS, 113, 367

Brightman, M., Harrison, F. A., Fürst, F., et al. 2018, NatAs, 2, 312

Carpano, S., Haberl, F., Maitra, C., & Vasilopoulos, G. 2018, MNRAS, 476, L45

Cenko, S. B., Barthelmy, S. D., D'Avanzo, P., et al. 2017, GCN, 21960

Coburn, W., Heindl, W. A., Rothschild, R. E., et al. 2002, ApJ, 580, 394

Deeter, J. E., Boynton, P. E., & Pravdo, S. H. 1981, ApJ, 247, 1003

Doroshenko, V., Tsygankov, S., & Santangelo, A. 2018, A&A, 613, A19

Epili, P., Naik, S., Jaisawal, G. K., & Gupta, S. 2017, MNRAS, 472, 3455

Ferrigno, C., Farinelli, R., Bozzo, E., et al. 2013, A&A, 553, A103

Finger, M. H., Bildsten, L., Chakrabarty, D., et al. 1999, ApJ, 517, 449

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd ed.; Cambridge: Cambridge Univ. Press)

Fürst, F., Pottschmidt, K., Wilms, J., et al. 2014, ApJ, 780, 133

Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, ApJL, 831, L14

Ge, M., Zhang, S., Lu, F., et al. 2017, ATel, 10907

Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Hanuschik, R. W. 1989, Ap&SS, 161, 61

Hastings, W. K. 1970, Biometrika, 57, 97

Israel, G. L., Belfiore, A., Stella, L., et al. 2017a, Sci, 355, 817

Israel, G. L., Papitto, A., Esposito, P., et al. 2017b, MNRAS, 466, L48

Jahoda, K., Markwardt, C. B., Radeva, Y., et al. 2006, ApJS, 163, 401

Jahoda, K., Swank, J. H., Giles, A. B., et al. 1996, Proc. SPIE, 2808, 59

Jaisawal, G. K., & Naik, S. 2016, MNRAS, 461, L97

Jaisawal, G. K., Naik, S., & Chenevez, J. 2018, MNRAS, 474, 4432

Jenke, P., & Wilson-Hodge, C. A. 2017, ATel, 10812

Jenke, P., Wilson-Hodge, C. A., & Malacaria, C. 2018, ATel, 11280

Jenke, P. A., Finger, M. H., Wilson-Hodge, C. A., & Camero-Arranz, A. 2012, ApJ, 759, 124

Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303

Kennea, J. A., Lien, A. Y., Krimm, H. A., Cenko, S. B., & Siegel, M. H. 2017, ATel, 10809

Klochkov, D., Staubert, R., Santangelo, A., Rothschild, R. E., & Ferrigno, C. 2011, A&A, 532, A126

Klus, H., Ho, W. C. G., Coe, M. J., Corbet, R. H. D., & Townsend, L. J. 2014, MNRAS, 437, 3863

Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502

Koliopanos, F., & Vasilopoulos, G. 2018, A&A, 614, A23

Kouroubatzakis, K., Reig, P., & Andrews, J. A. Z. 2017, ATel, 10822

Krimm, H. A., Holland, S. T., Corbet, R. H. D., et al. 2013, ApJS, 209, 14

Lamb, F. K., Pethick, C. J., & Pines, D. 1973, ApJ, 184, 271

Laplace, E., Mihara, T., Moritani, Y., et al. 2017, A&A, 597, A124

Ludlam, R. M., Miller, J. M., Arzoumanian, Z., et al. 2018, ApJL, 858, L5

Malacaria, C., Klochkov, D., Santangelo, A., & Staubert, R. 2015, A&A, 581, A121

Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791

Monageng, I. M., McBride, V. A., Coe, M. J., Steele, I. A., & Reig, P. 2017, MNRAS, 464, 572

Müller, D., Klochkov, D., Caballero, I., & Santangelo, A. 2013, A&A, 552, A81

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A. 2017, MNRAS, 467, 1202

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015a, MNRAS, 447, 1847

Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015b, MNRAS, 454, 2539

Nesci, R. 2017, ATel, 10989

Nowak, M. A., Vaughan, B. A., Wilms, J., Dove, J. B., & Begelman, M. C. 1999, ApJ, 510, 874

Paul, B., & Naik, S. 2011, BASI, 39, 429

Postnov, K. A., Gornostaev, M. I., Klochkov, D., et al. 2015, MNRAS, 452, 1601

Prigozhin, G., Gendreau, K., Doty, J. P., et al. 2016, Proc. SPIE, 9905, 99051I

Reig, P., & Nespoli, E. 2013, A&A, 551, A1

Shapiro, S. L., & Lightman, A. P. 1976, ApJ, 204, 555

Stanek, K. Z., Kochanek, C. S., Thompson, T. A., et al. 2017, ATel, 10811

Stella, L., White, N. E., & Rosner, R. 1986, ApJ, 308, 669

Sugita, S., Negoro, H., Nakahira, S., et al. 2017a, ATel, 10813

Sugita, S., Negoro, H., Serino, M., et al. 2017b, ATel, 10803

Sutton, A. D., Roberts, T. P., Walton, D. J., Gladstone, J. C., & Scott, A. E. 2012, MNRAS, 423, 1154

Tauris, T. M., & van den Heuvel, E. P. J. 2006, in Compact Stellar X-ray Sources, ed. W. H. G. Lewin & M. van der Klis (Cambridge: Cambridge Univ. Press), 623

Tsygankov, S. S., Lutovinov, A. A., Churazov, E. M., & Sunyaev, R. A. 2007, AstL, 33, 368

Tsygankov, S. S., Lutovinov, A. A., Krivonos, R. A., et al. 2016, MNRAS, 457, 258

Tsygankov, S. S., Lutovinov, A. A., & Serber, A. V. 2010, MNRAS, 401, 1628

van den Eijnden, J., Degenaar, N., Russell, T., et al. 2017a, ATel, 10886

van den Eijnden, J., Degenaar, N., Russell, T., et al. 2017b, ATel, 10946

Wang, Y.-M. 1981, A&A, 102, 36

Walton, D. J., Bachetti, M., Fürst, F., et al. 2018, ApJL, 857, L3