Phase resolved spectrum of the Crab pulsar from NICER

Although rotation powered pulsars (RPPs) were first discovered at radio wavelengths, they emit mainly at the much higher X-ray and $\gamma$-ray energies. So observations at these high energies are required to understand their emission mechanism. Further, the radio emission of RPPs is believed to be due to a coherent radiation mechanism, which is relatively difficult to study. In contrast, the high energy emission is relatively simple – ultra relativistic charges traveling along magnetic field lines and emitting very narrow beams of radiation along the direction of their momentum. Therefore a realistic magnetic field line geometry (e.g. dipole field with swept back field lines close to the light cylinder) is sufficient to qualitatively reproduce the observed folded light curve (FLC) of RPPs. Such work gained speed after the launch of the Fermi Gamma-Ray Space Telescope in 2008 (Atwood et al., 2007); see also Harding (2016).

The earliest magnetic field geometry used for modeling was the static dipole. More realistic magnetic fields were explored as time progressed – the vacuum retarded dipole, the aligned force-free magnetosphere, in which the magnetic and rotation axes are parallel, and the oblique force-free magnetosphere, in which they are not parallel. The force-free magnetospheres do not have any particle acceleration by definition. Therefore dissipative magnetospheres were invoked, that had a finite conductivity and allowed for some acceleration. However, the dissipative models did not model self-consistently the source of particles that provide the conductivity. This led to modeling of pulsar magnetospheres using particle-in-cell codes, which are still in their infancy; see Harding 2016 and references therein.

The magnetospheric models were coupled with different sites of particle acceleration, known as “gaps”, in which space charge deficiency leads to residual electric fields and therefore to particle acceleration. The popular gaps are the polar cap gap (PC), the outer gap (OG), the slot gap (SG) and its associated two pole caustic (TPC) model, and the separatrix layer (SL) model (see Harding 2016 and references therein).

Most of the above work was done at $\gamma$-ray energies. While each model explained qualitatively the FLCs of some RPPs, none of them explained the FLCs of all Fermi RPPs (Bai & Spitkovsky, 2010; Romani & Watters, 2010; Kalapotharakos et al., 2014; Harding, 2016; Pierbattista et al., 2016). Further, quantitative comparison of the observed and modeled FLCs (in terms of $\chi^{2}$) was often quite poor (Harding, 2016). However, several of the above models were able to reproduce the typical double peaked FLC even from a single magnetic pole of the RPP. They were able to reproduce the sharp peaks in terms of caustics, in which photons from different parts of the emission region arrived at the same rotation phase simultaneously, due to a combination of light travel time and special relativistic aberration (Cheng et al., 1986). However, analysis of the phase resolved spectrum (PRS) has lagged behind the analysis of the FLC in most of these models, until recently (Cheng et al., 2000; Brambilla et al., 2014).

The most popular RPP for such work is the Crab pulsar. It is a very luminous source, so it has been observed right from radio to $\gamma$-ray energies. Further, it is probably the only RPP for which the FLCs at most energies are not only similar in appearance, but are also aligned in rotation phase (except for the radio precursor; Abdo et al. 2010). This implies that the geometry of the emission region is common for most energies in the Crab pulsar. Therefore modeling its FLC at X-ray energies may be as useful as modeling it at $\gamma$-ray energies. This work obtains the high phase resolution (large number of bins per rotation period) FLC and PRS of the Crab pulsar at soft X-ray energies ($1-10$ keV) using the Neutron star Interior Composition Explorer (NICER) satellite observatory (Arzoumanian et al., 2014; Gendreau et al., 2016), which is the main result of this work. Here the focus will be mainly on the formation of caustics or cusps in the FLC of the Crab pulsar at soft X-ray energies.

Cheng et al. (1986) and Romani & Yadigaroglu (1995) were one of the earliest theorists to demonstrate the formation of sharp peaks in the $\gamma$ ray FLC of RPPs. Using a static dipole magnetic field, Cheng et al. (1986) showed that sharp peaks can form in the FLC due to photons from a large portion of the outer gap (OG) arriving at a common rotation phase (see their Figure $11$). These so called “caustics” or “cusps” form when one takes into account the combined effect of the special relativistic aberration and light travel time on a photon. Romani & Yadigaroglu (1995) demonstrated the same using the more sophisticated vacuum retarded dipole (VRD) magnetic field. Clearly the sweep back of the magnetic field lines in a VRD did not affect the process of caustic formation. These developments enabled the explanation of the well separated double peaks in $\gamma$ ray FLCs in terms of emission from a single magnetic pole of the RPP.

Then the Fermi Gamma-Ray Space Telescope obtained high quality $\gamma$ ray FLCs of several RPPs. This enabled statistical testing of more refined theoretical models of the formation of FLC. Romani & Watters (2010) tested three types of magnetospheric structures and two locations of the emission zone. The three magnetospheric structures were static magnetic field, modified VRD which they label “atlas” field, and “pseudo force free” (PFF) magnetosphere, which contains charges, but no currents. The two locations were the two pole caustic (TPC) and the outer gap (OG). They concluded that statistically OG was preferable to TPC. However, their best-fit values of $\alpha$, the angle between rotation and magnetic axis, and $\zeta$, the angle between rotation axis and line of sight, differ significantly from independently known values for some RPPs. In the work of Romani & Watters (2010) caustics are caused by piling up in phase of emission from different regions (altitudes) of the magnetosphere.

Bai & Spitkovsky (2010) use only the three dimensional force free magnetosphere (FF), which by definition has no particle acceleration; and they explore the TPC and OG locations. The FF magnetosphere has a conducting plasma, and the charges and currents are solved self-consistently. They conclude that both the TPC and OG models are unsuitable, but favor a modified version of OG which they call the “separatrix layer” gap (SL). In this model caustics form by an entirely different process – that of radiation from the same magnetic field line piling up in phase. This is possible because outside the light cylinder the magnetic field lines are almost radial, in the so called “split monopolar” geometry. They have named this kind of caustic formation the “sky map stagnation” (SMS) effect. These caustics form in the outer magnetosphere. Each magnetic pole gives rise to one peak in the FLC and the rise of flux to the peak as a function of rotation phase is usually faster than its drop after the peak.

There are two important aspects of caustic formation. First, it depends sensitively upon the geometry of the magnetic field lines, since that determines the amount of aberration and light travel time (Bai & Spitkovsky, 2010). Second, as emphasized by Bai & Spitkovsky (2010), the magnetic field structure at the light cylinder determines the shape of the polar cap on the surface of the neutron star, which in turn determines the shape and extent of the region of emission of photons, which affects the formation of caustics.

Kalapotharakos et al. (2014) use a dissipative magnetosphere, using different values of conductivity $\sigma$; when $\sigma=0$, the situation is similar to the VRD magnetosphere, while $\sigma=\infty$ implies the FF magnetosphere. For small values of $\sigma$ they find that the emission is mainly from the inner magnetosphere and the peaks in the FLC are not sharp. For large values of $\sigma$ the emission is mainly from the outer magnetosphere and the peaks can be sharp. For very large values of $\sigma$ the emission is entirely from the outer magnetosphere and is very non-uniform, so complex FLCs are obtained. They favor a hybrid model that is FF inside the light cylinder and dissipative outside (FIDO). In terms of caustic formation, the FIDO model is essentially two pole, like in the SG model, but without emission at low altitudes or at off peak phases, like in the OG model.

The above is a brief summary of the mechanism of caustic formation in high energy FLCs of RPPs. The results of this work can address some questions regarding caustic formation. To achieve this, two aspects have to be kept in mind. First, NICER consists of $56$ co-aligned detectors, each of which is essentially a small X-ray telescope. This clever design avoids (for most sources) two of the difficulties commonly faced by X-ray observatories – pile-up of photons and detector deadtime. Correspondingly, this design requires the calibration of $56$ detectors in terms of their long time light curves (LTLCs) and their spectral properties. This has been done in this work for the $49$ useful detectors of NICER.

Second, to obtain FLC and PRS at high phase resolution, one needs highly accurate rotation frequency $\nu$ and its time derivative $\dot{\nu}$ as a function of epoch, since these parameters of the Crab pulsar change significantly with time and since the NICER observations of the Crab pulsar occurred over a duration of about $2.35$ years. The $\nu$ and $\dot{\nu}$ have been obtained self consistently from the NICER data itself and only that data has been used for which the error in the derived rotational phase of each photon is less than $1/2048$ of a period. This enables one to obtain the FLC and PRS of the Crab pulsar at a resolution of $1024$ bins per period, and certainly better than $512$ bins per period. The details of how these $\nu$ and $\dot{\nu}$ are obtained, the treatment of the Crab pulsar glitch of $2017$ November, and the consistency with the radio ephemeris are explained in Vivekanand (2020). Typically the accuracy of $\nu$ obtained by NICER is $0.0009$ $\mu$Hz, with a standard deviation of $0.0006$ $\mu$Hz. The typical accuracy of the $\dot{\nu}$ is $4\times 10^{-14}$ Hz/s, with a standard deviation of $3\times 10^{-14}$ Hz/s.

Section $2$ discusses the observations and analysis of this work. The next three sections describe the calibration of the $49$ useful detectors of NICER. Sections $6$ and $7$ presents the spectrum for the Crab pulsar. The last section discusses the implications of this work.

NICER consists of $56$ co aligned detectors (Arzoumanian et al., 2014; Gendreau et al., 2016; Prigozhin et al., 2016). Its operating range is $0.2-12$ keV with a peak collecting area of $1900$ cm² at $1.5$ keV and below. Time resolution is $0.1$ microsecond ($\mu$s). At the time of doing this work, NICER had observed the Crab pulsar on $72$ different days, starting from $2017$ Aug $5$ to $2020$ Apr $27$, which is a duration of $2.73$ years, with multiple observations on $8$ of those days, resulting in $80$ observation identity numbers (ObsID). Ten of these ObsIDs had live times of less than $100$ sec, so they were not analyzed. Four of them were ignored for reasons specified in Vivekanand (2020), which gives details of the earlier observations spanning $1.72$ years and their analysis. The analysis of this work is slightly different because of the enhanced sensitivity requirement for spectral work.

The data were analyzed using NICER version $2020$-$04$-$23$_V$007$a software included in the HEAsoft distribution 6.27.2; the calibration setup was CALDB version XTI($20200722$). The analysis for each ObsID begins with the pipeline tool nicerl2 with the following parameters: angle between pointing and earth limb ELV $=30^{\circ}$, and angle between pointing and bright earth BR_EARTH $=40^{\circ}$ (Stevens et al., 2019; Bogdanov et al., 2019; Guillot et al., 2019; Miller et al., 2019; Trakhtenbrot et al., 2019); magnetic cut off rigidity COR_RANGE $\geq 4.0$ (Bogdanov et al., 2019; Guillot et al., 2019; Malacaria et al., 2019; Riley et al., 2019; Trakhtenbrot et al., 2019; van den Eijnden et al., 2020); and undershoot count rate UNDERONLY_RANGE $=0$-$50$; the rest of the parameters have default values (Stevens et al., 2019; Bogdanov et al., 2019). Next, the tool nimaketime is used to create good time intervals (GTI), using additionally the following two criterion: angle between pointing and sun SUN_ANGLE $>90^{\circ}$ (Bogdanov et al., 2019); and observations must be done when NICER in not in sun light, SUNSHINE $=0$ (Hare et al., 2020).

Next, the tool fselect is used to exclude data of detectors $14$, $34$ and $54$, whose data is known to be problematic (Bogdanov et al., 2019; Ray et al., 2019). Coupled with the detectors $11$, $20$, $22$ and $60$ that are permanently switched off at NICER, this leaves $49$ out of the original $56$ detectors for our analysis. Then, the original GTIs of the event file are replaced by the GTIs obtained above (if required, the FITS file extension name EXTNAME has to be changed back to GTI_FILT). Then fselect is used to filter out photon events outside the GTIs, using the gtifilter option. Next the LTLC is obtained using the tool extractor with a bin size of $2.5$ sec. This is used to exclude durations of observations in which the mean count rate is unusually low or high; this is done by editing the GTI extension of the event file and filtering using fselect with the gtifilter option.

Then, fselect is used to extract the data of each individual detector and extractor is used to obtain its LTLC. The tool fstatistic is used on the LTLC to obtain basic statistics of the count rate of each detector – the mean value, its standard deviation, the minimum and maximum values. Unusually low or high values indicate flaring, satellite not exactly on source, etc.; these durations are also excluded.

High background counts are identified by extracting a LTLC of bin size $8$ sec in the energy range $12-15$ keV and excluding those durations in which the count rate is statistically higher than $1.0$ (Bult et al., 2018). Details are given in Appendix A. NICER’s background is believed to be typically less than $1$ count/s/keV. It is obtained by processing data from several RXTE background regions. Background region $5$ of RXTE has an approximately power-law spectrum that is $<1.0$ counts/s/keV at $0.4$ keV and $<0.1$ counts/s/keV above $1.4$ keV (Keek et al., 2018). In the top panel of their Figure $2$, the background decreases to $\approx 0.01$ counts/s/keV at $10$ keV, while the flux of the Crab nebula plus pulsar is $3.7$ counts/s/keV at that energy, which can be derived from Figure 4. Thus the background is $\approx 370$ times weaker than the Crab nebula plus pulsar at $\approx 10$ keV.

High particle background due to geomagnetic storms are accounted for by the space weather parameter Kp (Ludlam et al., 2019; Miller et al., 2019), which is required to be below the value $5$. This is done by first converting the START and STOP times in the event files from the TT timescale to UTC dates and times and then checking the Kp values during the observation in the corresponding yearly WDC files¹¹1ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/wdc/yearly/. Details are given in Appendix B.

Finally, the event epochs that remain after the above stringent filtering are referenced to the solar system barycenter using the barycorr tool, using the JPL ephemeris DE430, with the position of the Crab pulsar at that epoch as input. The $\nu$ and $\dot{\nu}$, and their reference epoch for each ObsID are then used to compute the rotation phase of each photon, which is written into the event file. Vivekanand (2020) describes most of the above analysis in a different context. Later sections of this work describe the specific analysis required for the corresponding sections.

The filtering above reduced the total live time by $\approx 60$%, to $\approx 139$ ksec, leaving only $29$ ObsIDs with sufficient data to analyze, that also had their $\nu$ and $\dot{\nu}$ estimated by Vivekanand (2020), who does not implement this stringent filtering, who thus had many more photons for obtaining $\nu$ and $\dot{\nu}$. For those groups of ObsIDs that had sufficient number of photons in spite of the stringent filtering, it was possible to estimate $\nu$ and $\dot{\nu}$, which were consistent with those derived by Vivekanand (2020).

Figure 1 shows the LTLC of the Crab nebula plus pulsar using the $49$ useful detectors of NICER, in the energy range $0.2-12$ keV, for the $2.35$ years of observations used in this work, with a bin size of $10$ sec. This was obtained using the tool extractor. The mean count rate is $10384$ per second with a standard deviation of $75$, which is $0.72$% of the mean. The maximum and minimum values in Figure 1 are $10610.0$ and $10173.5$ counts per second, which implies that the peak-to-peak variation of the LTLC is about $\approx(10610.0-10173.5)/10384\approx 4.2$%. This is consistent with the recent measurement of the Crab nebula plus pulsar’s LTLC by the NuSTAR and Swift/BAT missions²²2(https://iachec.org/wp-content/presentations/2019/sessionI_-Madsen.pdf). This is also consistent with the findings in the recent past, that the reported decrease in the Crab nebula’s X-ray flux by $\approx 7$% over $\approx 2$ years starting from MJD $54690$ (mid $2008$), as reported by Wilson-Hodge et al. (2011), is no longer observed.³³3https://ntrs.nasa.gov/api/citations/20120015012/downloads/-20120015012.pdf⁴⁴4http://www.iucaa.in/iachec/talkmaterials//66/IACHEC-02-29-2016%20-%20Case.pptx In fact Wilson-Hodge et al. (2011) themselves stated that they cannot predict if the flux decline will continue in the future. This result is also confirmed in Fig. $1$ of Kouzu et al. (2013).

Figure 2 shows the statistics of the $49$ individual LTLCs. The dots with error bars show the mean count rate for each detector over the $2.35$ year duration. The crosses and the pluses represent the maximum and minimum values, respectively. All except detector $00$ show statistically similar values. Detector $00$ has a significantly lower mean count rate and a standard deviation slightly higher than the others. For example, the mean count rates of detectors $00$ and $01$ are $187.78$ and $224.02$, respectively, while their standard deviations are $8.15$ and $5.29$. The minimum values in Figure 2 are typically $-3.8(3)$ standard deviations below the mean, the error in the last digit given in the parenthesis; the maximum values are typically $+3.8(3)$ standard deviations above the mean. Now, a Gaussian probability distribution function has $0.000145$ of its area beyond the value $\pm 3.8$ standard deviations of the variable; this is obtained using the error function. This implies that half of that area, or $0.000073$, lies above and below $\pm 3.8$ standard deviations, respectively. This is consistent with the fact that for our sample size of $13454$ bins for each detector, each of $10$ s duration, one can expect one sample beyond $3.8$ standard deviations. Note that the live time for each detector in Figure 2 is $13454\times 10=134.54$ ksec and not the actual $\approx 139$ ksec because the extractor tool ignores bins of small fractional live time.

Figure 3 compares the LTLCs of the $49$ individual detectors with that of all detectors combined (shown in Figure 1), after normalizing each LTLC with its mean count rate. The $\chi^{2}$ per degree of freedom is below $1.5$ for $35$ detectors, between $1.5$ and $2$ for $7$ detectors and between $2$ and $2.25$ for $6$ detectors. The $\chi^{2}$ is derived using the formula $(y2-y1)^{2}/(y2+y1)$, where $y1$ and $y2$ are the counts (not count rate) in the corresponding time bins of the two LTLCs being compared. Thus, only Poisson statistical fluctuations are being accounted for in Figure 3. In such a complicated system such as NICER, it is reasonable to expect other kinds of fluctuations as well, not to mention minute source fluctuations also. Therefore, it is concluded that the LTLCs of all detectors are broadly similar.

One is justified in normalizing each LTLC with its mean count rate, because the variations of the latter in Figure 2 represent not merely statistical fluctuations, but something more like, say, differences in effective area. The error bar on the mean counts in that figure is the standard deviation, and not the error on the mean; this was done to better compare the mean value with the maximum and minimum values. For 49 detectors there are $49\times 48/2=1176$ distinct pairs, of which only $10$, $24$ and $43$ pairs have count rates differing by less than $3\sigma$, $5\sigma$ and $10\sigma$, respectively. So the mean count rate of each detector is statistically distinct from that of another detector, and most probably represents the difference in their effective areas.

Figure 4 shows the combined raw spectrum (uncalibrated) of the $49$ detectors, consisting of $\approx 139$ ksec of data from $29$ ObsIDs, from energy $0.2-12$ keV. The sharp dip in spectrum at around $2.2$ keV is due to the Gold (Au) edge⁵⁵5https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/nicer/docs/xti/-NICER-Cal-Summit-ARF-2019.pdf; the Gold is coated on the reflecting surface of the X-ray concentrators (XRC).

Figure 5 shows the comparison between Figure 4 and the corresponding spectrum of each detector, in terms of the $\chi^{2}$ per degree of freedom obtained over $1100$ channels in the energy range $1-12$ keV, obtained using the same formula as in Figure 3. The lower energy range was ignored because here the $\chi^{2}$ was very high, mostly due to the Carbon, Nitrogen and Oxygen edges at $\approx 0.3-0.5$ keV, and also due to the fact that our analysis is done with data beyond $1$ keV. $41$ of the $49$ $\chi^{2}$ values are below $2.0$; the extreme value of $3.9$ belongs to detector $61$. Once again, one expects non-Poisson fluctuations to exist in the data. The spread in $\chi^{2}$ is much larger than their formal statistical errors. However most of the $\chi^{2}$ gets contribution from the $2.2$ keV Gold edge. Therefore, one does not expect better agreement than this. It is therefore concluded that the raw spectra of the $49$ detectors are statistically similar.

Figure 5 is not merely Figure 2 plotted in a different way; each is sensitive to different physical parameters. The latter is sensitive to time variations of effective area, deadtime, etc., while the former is sensitive to optical loading of the detectors, calibration of the gain of the receivers that assigns an energy channel to each photon event, etc.

Figure 6 explores whether the spectrum of each detector changes as a function of time. It displays the $\chi^{2}$ between data earlier to and later than the date $2018$ Nov $10$, which roughly divides the data into two halves. The $\chi^{2}$ is estimated using the same formula as in Figure 3, after scaling the two spectra with their respective live times, because the cut-off date does not divide the data into exactly two equal halves. Just as in the previous figure, the $\chi^{2}$ values in Figure 6 are reasonable. It is therefore concluded that the raw spectrum of each detector does not change as a function of time during the $2.35$ year of observations.

In both Figure 5 and Figure 6, most of the estimated $\chi^{2}$ are larger than their formal errors. An extreme position to take would be to conclude that no detector behaves like any other. However, the acceptance or rejection of a hypothesis based on $\chi^{2}$ depends upon the confidence level, which can be subjective in some contexts. In complicated systems such as NICER one expects large $\chi^{2}$ due to unknown or immeasurable factors. In principle, one can chose for further analysis only those detectors whose $\chi^{2}$ is less than, say, 2.0, and live with the corresponding reduced sensitivity. Here statistics can not help one – only experience, intuition and the context will. To cite an example, in the good old days radio astronomers considered the detection of a point source as positive if it was detected at a signal to noise ratio of $5$. One can ask why $5$, and not $6$ or $4$? The answer is that this number was arrived at through common experience and not through pure statistics.

Wilson-Hodge et al. (2018) describe the procedure for estimating the dead time for NICER. There are two sources of deadtime. First, after reception of an event, the detector is unable to process further events for some time, during which it processes the previous event. This duration is recorded in the DEADTIME column of the event files. Second, if the source is bright, or the background is flaring, data buffers may saturate in the Measurement/Power Units (MPUs), which may stall recording of further events until the MPUs are reset. This causes loss of GTI, or equivalently, deadtime on the source. In this work the former will be called detector deadtime while the latter will be called GTI deadtime. The detector deadtime is contributed by all events, and not just the X-ray events. In Vivekanand (2020) only the detector deadtime was taken into account and that too only for X-ray events.

The procedure of Wilson-Hodge et al. (2018) was repeated almost exactly for four ObsIDs – $1013010113$ and $1013010147$, having the shortest ($120$ s) and longest ($23364$ s) live times respectively, and $1011010201$ and $1013010122$ having intermediate live times ($1841$ and $3850$ s respectively). While Wilson-Hodge et al. (2018) used $100$ phase bins per period, and $1000$ phase bins for estimating GTI deadtime, the corresponding numbers here are $256$ and $1000$. The results are similar for all four ObsIds, so the detailed results for ObsID $1013010147$ alone will be presented here.

The top panel of Figure 7 shows the mean X-ray counts per detector as a function of phase for ObsID $1013010147$ (i.e., the FLC). This has been obtained using cleaned events (Wilson-Hodge et al., 2018). The middle panel shows the mean GTI per phase bin after subtracting the expected mean value of $91.263$ s, which is obtained by dividing the live time by the number of phase bins. This has been obtained from the GTIs of the individual MPUs. The last panel shows the mean deadtime per detector, obtained using all events (Wilson-Hodge et al., 2018). In each panel of Figure 7, the results are first obtained per MPU and then averaged.

Two points are noteworthy in Figure 7. First, in the middle panel, the difference in GTI per phase bin, between the observed and the expected values, has a mean value of $-0.1$ ms, and an rms of $0.8$ ms; the peak to peak variation is $4.25$ ms. The mean value of GTI differs from the expected value by $-0.0001/91.263\approx-0.0001$%, while the peak to peak variation corresponds to $0.00425/91.263\approx 0.005$% variation. This implies that the GTI per phase bin is almost constant across the FLC and that its observed value is very close to the expected value; and so there is no significant loss of GTI due to MPU buffer overflows. This is further supported by the fact that, had there been GTI losses, it would have been anti-correlated with the FLC in the top panel which appears not to be the case.

Second, the detector deadtime fraction in the last panel of Figure 7 in the off-pulse region (phase $<0.2$) is $0.50$%, while at the pulse peak it is $0.77$%. The corresponding numbers in Figure 1 of Vivekanand (2020) are $0.46$% and $0.73$%. Thus the detector deadtime fraction obtained using all events is only $\approx 5-9$% larger than that obtained using only X-ray events. However, these numbers are approximate. The correct method of estimating the deadtime fraction or count rate is to first estimate them per MPU, and then average (Wilson-Hodge et al., 2018), since the GTI per MPU as well as the number of detectors per MPU could differ across the MPUs.

This is done in Figure 8. The mean off pulse deadtime fraction is $0.497(1)$%, while that at the pulse peak is $0.772(7)$%; these numbers are very similar to the approximate values above.

The results of Figures $7$ and $8$ hold for data of the other three ObsIDs also, except for lower statistical significance on account of lower live times.

Therefore, if the GTI deadtime is non existent for the Crab pulsar, the difference between detector deadtime obtained using all events and that obtained using only X-ray events, must depend upon the difference in the number of these two event types. This is demonstrated in Table 1, which contains event statistics from the filter file of that ObsID; for example the filter file for ObsID $1013010113$ would be “ni$1013010113$.mkf” found in the directory AUXIL. The filter file contains the numbers of various kinds of events registered by NICER, logged once every second in time. In particular, it contains two columns labeled TOT_ALL_CNTS and TOT_XRAY_CNTS. The former contains the number of all types of events, while the latter contains the number of what are deemed to be X-ray events, but are not guaranteed to be so⁶⁶6https://heasarc.gsfc.nasa.gov/docs/nicer/mission_guide/. The second and third columns of Table 1 contain the sum of TOT_ALL_CNTS and TOT_XRAY_CNTS, respectively, for the duration within the GTI of the corresponding event file, for $49$ useful detectors. For example, the number of TOT_ALL_CNTS and TOT_XRAY_CNTS for ObsID $1013010113$ registered within the GTI are actually $1114762$ and $1094359$, respectively, for $52$ detectors; after deleting the contribution of detectors $14$, $34$ and $54$, found in the columns MPU_ALL_COUNT and MPU_XRAY_COUNT in the filter file, the numbers become $1049364$ and $1030023$, respectively.

The fourth column of Table 1 gives the ratio of the number of events of column $3$ that were eventually classified as X-ray events, obtained in the keyword NAXIS2 in the corresponding event file, and TOT_ALL_CNTS. Table 1 has been made for all $29$ ObsIDs used in this work, and the mean value of the above ratio is $0.922(5)$. This demonstrates that the difference between actual X-ray counts and all counts for the Crab pulsar data of NICER is typically about $7.8(5)$%, which is a small number.

There are two conclusions so far in this section. First, The difference between the number of all events and the number of X-ray events for the NICER data of the Crab pulsar is quite small ($\approx 8$%). This is not surprising, since the event rate for the Crab pulsar is dominated by events from the source; for much less luminous sources, the event rate would be dominated by particle events, etc. Second, the GTI deadtime is negligible for the Crab data of NICER. This is also not surprising, since source count rates of about $18846$ per second (for 49 detectors) begin to show such effects⁷⁷7https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/-nicer_analysis_tips.html#Deadtime, and the Crab pulsar’s peak count rate is only about ${{\approx 17000}}$ counts per second. Further, MPU buffer overflow depends upon the product of the count rate and the duration for which it is maintained;

Figure 9 compares the detector deadtime fraction of each of the $49$ detectors with the combined detector deadtime fraction of all $49$ detectors, across $256$ phase bins in a period. It estimates the $256$ differences of the two curves and computes their mean and standard deviation. Ideally, the mean should be zero for each detector and the standard deviation should reflect the average measurement error of the detector deadtime fraction across a period. Figure 9 shows the ratio of mean and standard deviation for each detector. The mean value of this ratio in Figure 9 is $-1.6\pm 2.0$. If one ignores the extreme value of $5.6$ of detector $62$, the results are $-1.7\pm 1.7$. Thus the mean value of this ratio is consistent with being zero.

In Figure 9, the relatively large values of standard deviations (not standard deviation of the ratio mean/rms) for each detector are due to two factors – the reduced number of photons in each phase bin for a single detector and, to a lesser extent, the intrinsic distribution of detector deadtimes in the column DEADTIME in the event files. It should be recalled that the quantity being estimated here (detector deadtime fraction) is very small, of the order of $0.5$% to $0.8$%; the difference of two such quantities is typically of the order of $\approx 0.05$% to $0.005$%. Now, the mean number of events in a phase bin for a single detector is $1570142182/49/256\approx 125171$; the Poisson fluctuation in this quantity is of the order of $\sqrt{(}125171)/125171\approx 0.28$%; over $256$ phase bins the accuracy of the estimated mean (due to the Poisson process alone) would be $0.28/16\approx 0.02$%, which is comparable to the quantity being estimated. Similarly, the mean detector deadtime for ObsID $1013010147$ is $23.1\pm 4.2$ $\mu$s (Vivekanand, 2020). The fractional error is $4.2/23.1\approx 18.2$%. The mean detector deadtime estimated using $125171$ events will have an accuracy of $18.2/\sqrt{(}125171)\approx 0.05$%, which is also comparable to the quantity being estimated.

The time evolution of detector deadtime is studied by estimating the deadtimes of the first and second halves of the data, for all $49$ detectors combined, with the cutoff date of $2018$ Nov $10$ that was used in Figure 6. The quantity plotted in Figure 10 is $(z2-z1)/(z2+z1)*2.0$ as a function of phase, where $z1$ and $z2$ are the detector deadtime fractions in the first and second halves of the data, respectively. The mean value of the fractional difference is $-0.016(1)$; i.e., the earlier and later detector deadtime fractions differ by about $1.6(1)$%. It is therefore concluded that there is no evolution of detector deadtimes as a function of epoch, during the $2.35$ years of observations studied here.

Finally, Figure 11 investigates the dependence of detector deadtime upon photon energy. The combined data of all $49$ detectors was divided roughly into two equal parts – one consisting of photons with energy between $1.0-1.62$ keV, and the other between $1.62-10.0$ keV. The FLC and deadtime fraction for each energy band is shown in the top and bottom panels of Figure 11, respectively. The spectral evolution of the FLC with energy is evident in the top panel of Figure 11. The deadtime fraction curves mirror this evolution in the bottom panel, so it is not possible to directly compare the two curves. However, if the detector deadtime in a phase bin depends only upon the photon count rate in that bin and not upon the photon’s energy, then the deadtime fraction would have the same linear dependence upon the count rate for the two energy bands. So a linear fit was done to the deadtime fraction against the count rate in each energy band (see Vivekanand (2020) for justification of linear dependence), and the results are given in Table 2.

The two slopes are almost the same, while the intercepts are negligible. This proves that the detector deadtime fraction depends only upon the count rate and not upon the photon energy.

The top panel of Figure 12 shows the calibrated phase averaged spectrum (PAS) of the Crab nebula plus pulsar using all $49$ useful detectors, for the $29$ ObsIDs with a live time of $\approx 139$ ksec. The spectrum was obtained using the extractor tool. It was analyzed using XSPEC⁸⁸8https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/manual.html version $12.11.0$, using nixti20170601_combined_v002.rmf as the main response file and nixtionaxis20170601_combined_v004.arf as the ancillary response file. These differ from the corresponding files found in CALDB version XTI($20200722$), viz. nixtiref20170601v002.rmf and nixtiaveonaxis20170601v004.arf, in that the latter two are relevant for observations using all $52$ detectors of NICER, while the former two are relevant for observations using only $49$ detectors; see the relevant analysis thread⁹⁹9https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/arf-rmf/ at the NICER site. Channels $230$ to $1000$ were used (energy range $2.3-10$ keV). No background spectrum was used since that becomes important for the Crab pulsar only at much higher energies (Madsen et al., 2015), as mentioned in section $2$. Ludlam et al. (2018) state that the background fields have count rates of $0.7-2.3$ counts/s, while the Crab pulsar’s average count rate is $\approx 11020$ counts/s, for $52$ detectors. The default photo ionization cross-sections and solar abundances were used (vern and angr). An absorbed power law model was fit to the data of Figure 12; the parameters of the model are the column density of absorbing Hydrogen $N_{H}$, the power law index $\Gamma$, and the power law normalization $A$. The results are: $N_{H}=0.210(3)\times 10^{22}$ cm^-2, $\Gamma=2.050(1)$; and $A=8.82(1)$; the $\chi^{2}$ per degree of freedom is $1.99$, for $769$ degrees of freedom. The errors are obtained using the error function and represent $90$% confidence levels; at higher decimal places the lower and upper limits of the errors are asymmetric.

These results are consistent with the results of XMM-Newton in the $0.6-9$ keV band for the overall Crab region (Kirsch et al., 2006), in their Table $2$ and Figure $6$. Their values are $N_{H}=0.260(1)\times 10^{22}$ cm^-2 and $\Gamma=2.046(3)$, with $\chi^{2}=1.99$ per degree of freedom.

These results also compare well with the results of NuSTAR (Madsen et al., 2015) in their Figure $1$, in the energy range $3-78$ keV. They fixed the value of $N_{H}$ at $0.2\times 10^{22}$, since their data is not sensitive to this parameter, and use the wilms solar abundances. Their $\Gamma=2.0963(4)$, and they do not mention the $\chi^{2}$. Their $\Gamma$ is in the same numerical range as our result considering the fact that their energy range is much larger.

The bottom panel of Figure 12 shows the ratio of the data to the absorbed power law model curve. As found in Figure 1 of Madsen et al. (2015), the variations are of the order of $1$%, presumably indicating the accuracy of the spectral calibration of these instruments.

As mentioned by Madsen et al. (2015), the non piled up instruments covering the $1-100$ keV band agree on a power law index of $2.10(2)$ for the PAS of the Crab nebula plus pulsar, without any break or curvature in this band; see also Kirsch et al. (2005). Kouzu et al. (2013) also confirm this result using the HXD detector aboard Suzaku. Our result is consistent with this statement. Moreover, the more important result is the variation of $\Gamma$ across the rotation phase, and not its absolute value, since that is correlated with $N_{H}$ (see Figure $2$ of Massaro et al. 2000; see also XSPEC manual).

It is therefore concluded that the PAS of the Crab pulsar and nebula in Figure 12 is consistent with earlier studies.

Figure 13 is similar to Figure 12, except that lower energy limit is $1.0$ keV. The ratio of the data to the model in its lower panel has variations of the order of $1$%, the same as in the latter figure. The fit parameter values are: $N_{H}=0.2943(1)\times 10^{22}$ cm^-2, $\Gamma=2.0651(2)$; and $A=9.095(2)$. These are consistent with the values obtained above for Figure 12, but the errors are unreliably smaller, because the $\chi^{2}$ per degree of freedom is $44.4$, contributed mainly by the channels between energies $1.0-2.3$ keV. This appears to be due to the very high photon counts in the lower energy channels. This problem disappears when the PRS is computed in the next section, since the number of photons in each energy channel is now less by the number of phases in the period.

The PRS of the Crab pulsar has been studied over the last two and half decades by several workers at both the soft and the hard X-ray energies, using different X-ray observatories (RXTE, BeppoSAX, XMM-Newton, Chandra, NuSTAR and Insight-HXMT) (Pravdo et al., 1997; Massaro et al., 2000, 2006; Kirsch et al., 2006; Weisskopf et al., 2011; Ge et al., 2012; Madsen et al., 2015; Tuo et al., 2019). The highest phase resolution achieved so far has been typically $\approx 1/100$ in phase; even though BeppoSAX (Massaro et al., 2000, 2006) and RXTE (Ge et al., 2012) start with $300$ and $1000$ phase bins for the FLC, respectively, they have had to estimate the PRS over $\approx 100$ phase bins to maintain statistical significance. The following sub section presents the PRS of the Crab pulsar obtained by NICER over $128$ phase bins. The next sub section obtains the PRS over $1024$ bins, although the number of independent phase bins is probably $512$ or better.

The $29$ ObsIds used in this work already have their $\nu$ and $\dot{\nu}$ (each at its reference epoch) estimated by Vivekanand (2020), along with their formal errors. These were verified in those cases where the event files had sufficient number of photons in spite of the stringent filtering of this work. For each ObsID, using its start and stop times available in the corresponding event file, a phase error was calculated, that would result if the $\nu$ and $\dot{\nu}$ used were different by twice their formal errors. Only those ObsIDs were retained that had the phase error less than $1/2048$. This eliminated the ObsIDs $1011010201$ and $1013010119$, leaving $27$ ObsIDs for further analysis. Next, FLCs were formed for each of the $27$ ObsID in $1024$ phase bins. These were then cross correlated with the FLC of ObsID $1013010147$, after shifting it in phase so that the first peak of its FLC lies at phase $0.4$. The $26$ phase offsets thus obtained had a typical error of $\approx 10^{-4}$ or better. These phase offsets were then used along with the corresponding $\nu$ and $\dot{\nu}$ to compute the phase of each photon event in the $27$ ObsIDs (see Vivekanand (2020) for details).

This section will begin with the scientific motivation for obtaining high resolution FLC and PRS of the Crab pulsar at high energies. It was mentioned in section $1$ that these are required to understand the high energy emission mechanism of RPPs. It was also mentioned that the analysis of the PRS lagged far behind that of the FLC. One reason was probably the success obtained in theoretically producing the double peaked FLCs of several RPPs based only on the geometry of the magnetic field. The other reason was the relative difficulty in applying the basic radiative processes in a realistic magnetosphere of the RPPs, viz., curvature radiation (CR), inverse Compton scattering (ICS) and synchrotron radiation (SR) (Cheng et al., 1986). The availability of high resolution PRS, and its use in conjunction with high resolution FLC, necessitates the invocation of these basic radiative processes in addition to magnetic field line geometry and magnetospheric processes. Such a problem is very complicated. In this work, the focus is restricted to the formation and properties of caustics in FLCs using the main results of this paper (Figure 15).