A Comprehensive Library of X-ray Pulsars in the Small Magellanic Cloud: Time Evolution of their Luminosities and Spin Periods

We have collected and analyzed the complete archive of {\itshape XMM-Newton\} (116), {\itshape Chandra\} (151), and {\itshape RXTE\} (952) observations of the Small Magellanic Cloud (SMC), spanning 1997-2014. The resulting observational library provides a comprehensive view of the physical, temporal and statistical properties of the SMC pulsar population across the luminosity range of $L_X= 10^{31.2}$--$10^{38}$~erg~s$^{-1}$. From a sample of 67 pulsars we report $\sim$1654 individual pulsar detections, yielding $\sim$1260 pulse period measurements. Our pipeline generates a suite of products for each pulsar detection: spin period, flux, event list, high time-resolution light-curve, pulse-profile, periodogram, and spectrum. Combining all three satellites, we generated complete histories of the spin periods, pulse amplitudes, pulsed fractions and X-ray luminosities. Some pulsars show variations in pulse period due to the combination of orbital motion and accretion torques. Long-term spin-up/down trends are seen in 12/11 pulsars respectively, pointing to sustained transfer of mass and angular momentum to the neutron star on decadal timescales. Of the sample 30 pulsars have relatively very small spin period derivative and may be close to equilibrium spin. The distributions of pulse-detection and flux as functions of spin-period provide interesting findings: mapping boundaries of accretion-driven X-ray luminosity, and showing that fast pulsars ($P<$10 s) are rarely detected, which yet are more prone to giant outbursts. Accompanying this paper is an initial public release of the library so that it can be used by other researchers. We intend the library to be useful in driving improved models of neutron star magnetospheres and accretion physics.


INTRODUCTION
The purpose of this paper is to motivate and facilitate the community to investigate on a population-statistical basis, the processes of accretion and X-ray emission in High Mass X-ray Binary (HMXB) pulsars. To this end we have applied a uniform reduction/analysis pipeline to compile a comprehensive set of physical parameters from the large archive of observations of the Small Magellanic Cloud (SMC).
The SMC is a dwarf irregular galaxy at a distance of 62 kpc from the Milky Way (Graczyk et al. 2014;Scowcroft et al. 2016). It contains a large and active population of HMXBs (e.g. Galache et al. 2008;Townsend et al. 2011;Klus et al. 2014;Coe & Kirk 2015;Haberl & Sturm 2016;Christodoulou et al. 2016). Systems involving a Be type star account for 98% of the confirmed SMC HMXBs .
Be/X-ray binaries (Be-XBs) are stellar systems in which a neutron star (NS) accretes matter from the circumstellar disk of a massive early-type companion resulting in outbursts of high-energy radiation. During all such outbursts the X-ray flux is pulsed at the spin period of the neutron star, due to magnetic channeling of the accretion flow to the poles of the neutron star, hence these objects are transient X-ray pulsars. The triggering of outbursts is partly controlled by the orbital parameters of the binary, with many systems exhibiting strings of sub-Eddington outbursts spaced at the orbital period, occurring close to periastron, and lasting for about 1/10 of the orbit. Not all periastron passages lead to an outburst so other factors must be involved, and are the subject of current research. Some Be-XBs have undergone infrequent but much more luminous outbursts, exceeding the Eddington luminosity (L Edd ) and lasting for a complete orbit or longer (e.g., Martin et al. 2014). The cause of the latter "giant outbursts" is suspected to be due to a sudden increase in the mass-transfer rate from the companion, driven by tidally induced density perturbations in the circumstellar disk (Negueruela et al. 1998;Moritani et al. 2013). These evolving dynamic structures are observed in the Hα line profiles of Be-XBs (Reig et al. 2016). It is also possible that a radiative instability in the pulsar's accretion disk, similar to the situation leading to a nova eruption, plays a key role in triggering giant outbursts .
Proximity, compact size, and minimal foreground contamination have made the SMC an ideal laboratory to study HMXBs. The motivation of this paper is to advance the study of its population as a whole, to enable robust statistical analysis of the physical parameters of HMXBs, and in turn enable a new generation of theoretical models to improve understanding of accretion physics. The X-ray source population of the SMC is different than that found in the Milky Way, the Large Magellanic Cloud (LMC), and other Local Group galax-arXiv:1703.05196v1 [astro-ph.HE] 15 Mar 2017 ies (e.g., M31, M33). Many more HMXBs are known in the SMC than in the LMC and the Milky Way (Haberl et al. 2000;Yokogawa et al. 2003). Based on the masses of the Milky Way and the SMC, the SMC HMXBs are a factor of 50 more numerous than what one would expect (Walter et al. 2006;McBride et al. 2008). The SMC is experiencing an era of ongoing star formation and its large number of HMXBs is certainly related to the high star formation rate (SFR) (Grimm et al. 2003). The timescale for production of HMXBs in the SMC has been revealed by Antoniou et al. (2010) to peak in the 25-60 Myr range, making the pulsars in our sample of similar age. The above factors establish the SMC as a unique place to study this important branch of stellar evolution.
Building on the legacy of almost 2 decades of dedicated monitoring of the SMC, in this paper we study observations collected in surveys carried out with XMM-Newton (Haberl et al. 2008), Chandra (Antoniou et al. 2009;Laycock et al. 2010), and RXTE (Laycock et al. 2005;Galache et al. 2008). The basis of our analysis is the SMC X-ray pulsar (SXP) catalog of Coe & Kirk (2015). This catalog lists the coordinates, orbital period, eccentricity, measured spin period of the compact object, and characteristics of the companion Be star where known for currently identified SMC X-ray pulsars. We have mined the data from the archives of these three telescopes and we have generated a new comprehensive time-domain library of high-level data products for the SMC pulsars in our resulting sample.
This new archive incorporates all previous surveys of the SMC and adds all the latest observations that fall in the public domain up to the year 2014. For each known SMC pulsar, all of these archived data are combined together in order to produce a complete picture of the X-ray emission. The library also contains comprehensive (folded and unfolded) light curves at different energy bands, the variations of the luminosities, the pulsation amplitudes (count rate), the pulsed fractions, anḋ P information for each known pulsar, all of which portray the long-term behaviors of these objects.
The resulting products can be used toward making progress in the following areas: (a) map out the duty cycles of X-ray emission from HMXBs and delineate the various phases of accretion and quiescence; (b) provide a library of pulse profiles to confront geometric models (e.g. Yang et al. 2017) of pulsar emission and constrain their physical and geometrical parameters; (c) investigate statistical correlations between the physical properties of the compact objects, across the full parameter space; and (d) produce large-number statistics for the entire class of objects. This wealth of information will also be released for public use. Release 1.0 (coincident with the publication of this paper) comprises the catalog of measurements obtained from our pipeline processing. The scope and contents of the available database is illustrated by the plots and tables presented throughout this paper. Release 2.0 will include pulse profiles and photon event files.
In the following sections we describe the contents of the new library, how the information is grouped, and how it can be used. In § 2, we describe the observations from the three observatories and our processing of the raw data from each archive. In § 3, we describe the products and content of the library, illustrated with examples of individual pulsars, and present some statistical inferences from the combined data set of the 3 satellites and SMC pulsars. We conclude in § 4 with a discussion of our products and a summary of the paper.

OBSERVATIONS
The library is constructed from a large volume of archival XMM-Newton, Chandra and RXTE data, which is summarized in Table 1 and described in detail in the following subsections.
2.1. XMM-Newton Observations XMM-Newton was launched in December 1999 by the European Space Agency. It carries three identical grazing incidence X-ray telescopes, whose point spread function (PSF) is sufficient to spatially separate most of the individual bright X-ray sources of the SMC. Up until 2014 (the cutoff for this project), 116 XMM-Newton observations of the SMC were available in the archive, as listed in Table 1.
XMM-Newton has the largest effective area of any Xray satellite at energies below 2 keV, a record it will retain until the Neutron Star Interior Composition Ex-ploreR (NICER) launches; it also has a 30 cm optical/UV telescope (the Optical Monitor; Mason et al. 2001) allowing simultaneous X-ray and optical/UV coverage. XMM-Newton's X-ray telescopes each feed one of the European Photon Imaging Cameras (EPIC) (Strüder et al. 2001;Turner et al. 2001), which comprise the PN, and Metal Oxide Semi-conductor (MOS1, MOS2) instruments covering the energy range 0.1-15 keV with moderate spectral resolution.
From the XMM-Newton archive we acquired the EPIC PN data (which have a higher time resolution than the MOS data) for all publicly available SMC observations obtained from 2000 to 2014. In Fig. 1 the number of observations in which the different known SMC pulsars appeared in the field of view (FOV) and the number of positive detections of them by the PN camera are shown in green and red, respectively. The spin periods of these SMC pulsars range from 0.72 to 4693 s. Here positive detections means photon counts and flux are recorded in the XMM-Newton Science Archive (XSA) 6 . The Xray source detections are above the processing likelihood of 6 7 . In this histogram, we can see how many times each source was caught in quiescence (corresponding to non-detections) and we can identify the sources that are appear to be permanent (at least in the context of our 15 year study baseline) rather than transient emitters.
In our data analysis pipeline, we begin with the SXP catalog of Coe & Kirk (2015) which contains pulse periods and celestial coordinates of each known pulsar, determined from the existing body of publications which includes X-ray and optical counterpart identifications. Most of the SXP pulsars have positions known to subarcsecond accuracy. We take the XMM-Newton source catalog obtained from the XSA and search for EPIC detections (both MOS and PN) in a positional search radius for each known pulsar. The initial search uses a search radius of 15 about the known SXP coordinates. Every XSA catalog point-source within this radius is evaluated based on its positional offset and uncertainty, requiring the offset to be smaller than 3σ c for a positive identification with the SXP object, where σ c is the combined uncertainty computed via Equation 1, and where r of f is the offset between the known position and the detected point source position, r xmm is the XMM-Newton Science Archive (XSA) Right Ascension and Declination (RADEC) combined error (determined while fitting the detection), r sys is the XSA systematic error of the XMM-Newton fields, and r psr is the uncertainty of the known position of the pulsar. Our data reduction of these sources included the standard procedures of the XMM-Newton Science Analysis Software (SAS, version 1.2) from the XMM-Newton Science Operations Center (SOC). We retained the standard XMM-Newton event grades (patterns 0-4 for the PN camera), as these have the best energy calibration. The PN data were analyzed using the commands evselect and epiclccorr in SAS. For the EPIC detectors, (Strüder et al. 2001;Turner et al. 2001), the data for the pointlike sources were extracted from circular regions of radius 20 . Background levels were estimated and subtracted using annular regions defined by inner and outer radii of 30 and 60 centered on each source.
Source fluxes were directly obtained from the XSA using the known pulsar coordinates. If a source was not detected but it was in the FOV, then we recorded the upper limit to the flux that was calculated from the "Flux Limits from Images from XMM-Newton" (FLIX) 8 server.
It was necessary to extract XMM-Newton light curves from scratch at the native timing resolution of the detectors (0.0734 s for PN) because the XSA pipeline versions are generated using a binning scheme that precludes high-resolution timing analysis. We also extracted and archived the event list (Time, Energy) for each individual source. This product is the starting point for performing a range of advanced event-based analyses, e.g. FFT (Israel & Stella 1996), quantiles (Hong et al. 2005), energy dependent pulse profiles and 2D phase-energy-intensity histograms (Schönherr et al. 2014;Hong et al. 2016 The times of arrival of the photon events were shifted from the local satellite frame to the barycenter of the solar system using the task barycen in the SAS software. Additional improvements using the cleaned event file and the 'gti' (good time interval) task tabgtigen were also made in order to exclude the times when the background was very high. The rate threshold of the 'gti' fits file applied to the event file was set at 0.6 counts/s. A pulsation search was performed following the procedure described in Section 2.4.1 to search for periodicities in the broad, soft (0.2-2 keV), and hard (2-12 keV) energy bands.

Chandra Observations
The Chanda X-ray Observatory was launched in July 1999, carrying the High Resolution Mirror Assembly (HRMA), the highest angular-resolution X-ray optics ever flown in space. The HRMA focusses X-ray light on either the Advanced CCD Imaging Spectrometer (ACIS) or the High Resolution Camera (HRC). For this project we used ACIS data, which provide sub-arcsecond source positions, photon energies, and event timing to 3.2 s resolution in standard operation.
The known X-ray pulsars in the SMC (Coe & Kirk 2015) were searched for in the first fifteen years of Chandra observations obtained from the Chandra Data Archive. We applied the latest Chandra calibration files to each image and we reduced the data with the Chandra Interactive Analysis of Observations software package (CIAO, version 4.5) (Fruscione et al. 2006). We created X-ray images within the 0.3-7.0 keV energy band in which Chandra is best calibrated and most sensitive. First, the image files were obtained by executing fluximage on the reprocessed event files. We then used mkpsfmap to generate a corresponding image whose amplitude is the PSF size in terms of the region that encloses 95% of the counts at 1.5 keV. The task evalpos was used to look up the PSF size at each source position in this image, and dmmakereg was used for generating the appropriate sized circular source extraction and annular background extraction regions.
After verifying the accuracy of the Chandra astrometry for the fields of interest, we placed our extraction regions at the coordinates of each pulsar to extract the source and background events. Finally, source fluxes and light curves were extracted with the CIAO tools srcflux and dmextract, respectively. The tool srcflux was used with a power-law spectral model with photon index 1.5 and an absorption column 5 × 10 21 cm −2 , representative of the  line-of-sight to the SMC.
srcflux was also used to determine whether the pulsar was detected, in which case a flux measurement is reported, or not, in which case an upper limit is reported instead. The source flux was extracted within the region that encloses 95% of the counts at 1.5 keV. This region was determined for each detection seperately using the PSF map. The srcflux derived net photon counts 1 at 90% (1.65 σ) confidence level were recorded as detections. If the net counts are 0, then the upper limit is reported. For srcflux positive detections, a source event file, light curve, and pulse-height spectrum were extracted, and saved in the library along with the appropriate response files. Light curves were extracted at 3.2 s resolution, set by the read-mode most commonly encountered in the archival dataset. A pulsation search was performed following the procedure described in Section 2.4.1.
The number of times each source was in a Chandra FOV and the number of times each source was detected by Chandra are shown in Fig. 2. As in the analogous Fig. 1 for the XMM-Newton data, this plot serves to illustrate the frequency of each source being caught in quiescence (non-detections) and allows us to identify the "permanent" rather than transient emitters.

RXTE Observations
The Rossi X-ray Timing Explorer (RXTE ) launched in December 1995 carrying the Proportional Counter Array (PCA), All-Sky Monitor (ASM), and High Energy X-ray Timing Experiment. The PCA consisted of 4 Xenonfilled Proportional Counter Units (PCUs), and provided individual event timing at µs resolution and moderate energy resolution over the range 2-60 keV. Starting in 1997 the SMC was observed approximately weekly with the PCA for some 16 years, accumulating ∼1000 observations (Table 1) before RXTE was deactivated in 2012 January. Data in Good Xenon mode were extracted in the 3-10 keV energy range (maximizing S/N) using the FTOOLS suite. In our subsequent analysis, we used the output of the IDL pipeline products (PUlsar Monitoring Algorithm or PUMA) up to year 2012 (Galache et al. 2008). To summarize the data reduction: the data were first cleaned using standard FTOOLS scripts 9 , then the FTOOL maketime was used to generate the GTI files. We extract the light curves with 0.01 s binning 10 and we applied background subtraction and barycentric correction. Finally, the count rates were normalized to account for the varying number of active Proportional Counter Units. Further details can be found in Galache 2006;Galache et al. 2008;Townsend et al. 2011Townsend et al. , 2013Klus et al. 2014.
The number of RXTE observations with collimator response > 0.2 of each known SMC pulsar is shown in Fig. 3. RXTE PCA was not an imaging detector, and this figure shows only the number of observations for which the coordinates of each pulsar were in the 2 • Full Width at Zero Intensity (FWZI) FOV. Following the approach established in prior works based on the Lomb-Scargle (LS) periodogram (Laycock et al. 2005(Laycock et al. , 2010Galache et al. 2008), we searched for pulsations in the light curves of each pulsar from each satellite detection. The significance, s, of each periodicity was calculated from the number M of independent frequencies and the LS power P X according to Press et al. (1992), (2) The error in spin period was calculated from the standard deviation of the frequency (Horne & Baliunas 1986), where σ 2 N is the variance of the light curve, N is the number of data points, T is the total length of the data, and A s is the amplitude of the signal.
Each periodogram was automatically scanned to look for the expected fundamental and harmonics of the pulsar. For the imaging instruments the specific pulsar is known, and so the search is very targeted. For RXTE the periods of all known pulsars in the FOV of the PCA collimator are searched, and we used the fundamental harmonics of [0.8 − 1.2] × P with the s 99% and all the collimator responses in the later analysis. Confidence levels are assigned following the prescription of Press et al. (1992) computed for a search over a 10% tolerance range on the expected period. We treat significance 99% as a valid detection and do not adjust the threshold to account for the number of observations. This ensures a uniform criterion that does not change if the definition of the sample were to change. When a detection was made we obtained the spin period, then folded the light curve and obtained the pulsation amplitude and the pulsed fraction (PF). The pulsed fractions of the light curves were calculated by integrating over the pulse profile according to the prescription of Bildsten et al. (1997). Here where n bin is the number of bins for each folded light curve, f j,mean is the mean photon count rate in each bin, f min is the minimum of f j,mean , and f i is the photon count rate of the un-binned light curves. Note that for RXTE we did not compute pulsed fractions since the PCA is a non-imaging detector and multiple sources are always in the FOV, so the un-pulsed component cannot be reliably measured.

Cross Calibration
In order to combine the data for the three satellites together it is necessary to account for differences in sensitivity and energy range. In principle the response functions of each instrument are well characterized and good cross calibration for flux and luminosity is possible for the imaging instruments. Since RXTE is not an imaging instrument, and has very different properties, we took an empirical approach to scale the pulse  amplitudes. For the purposes of plotting our results in this paper we performed a cross-calibration to normalize the luminosities and pulse amplitudes (count rates) to be on the same scale as the XMM-Newton PN. In the case of Chandra this was accomplished using the tool PIMMS 11 (Portable, Interactive, Multi-Mission Simulator) to convert the count rate from the ACIS-I detector (with no grating) into XMM-Newton PN count rate (with no grating and medium filter). The input energy range was 0.3-7 keV and the output energy range was 0.2-12 keV. We used power-law models with absorption 5 × 10 21 cm −2 (calculated using NASA's HEASARC tool 12 for the SMC sightline) and photon indices of 1.5. With a 1 count/second Chandra detection, the predicted XMM-Newton count rate is about 4 counts/s. Thus, the pulsation amplitudes from Chandra were scaled by a factor of 4.
For the RXTE data we do not directly measure absolute flux due to the likelihood of source confusion. Many observations feature 2 or more pulsars active simultaneously (commonly 1-3, but in one case 7). Instead we followed an empirical approach to estimate total flux from pulse amplitude (which we measure directly): we used PIMMS to calculate the absorbed fluxes from the measured pulse amplitudes, assuming a fixed pulsed fraction. We set both the input and the output energy range of the PCA detector to 3-10 keV. From the PIMMS power-law model described above, the predicted flux for 1 count/PCU/s was F = (9.23 −0.12 +0.10 )×10 −12 erg cm −2 s −1 . A proxy for the X-ray luminosity was then calculated according to where A is the pulsation amplitude in the unit of counts/s, P F is the (unknown) pulsed fraction which can vary between 0.1 and 0.5 , and r = 62 kpc is the distance to the SMC. In the last step of eq. (5), we chose to set P F = 0.4 so that our X-ray luminosities are comparable with the corresponding values obtained by Klus et al. (2014) who used average count rates (over many observations) instead of individual A values to calculate fluxes.

OVERALL PROPERTIES OF THE SMC LIBRARY
Our library includes data products from the XMM-Newton, Chandra and RXTE satellites processed using our data reduction and analysis pipelines, together with parameters computed from the combined products. The number of times that each pulsar was observed by each satellite is reported in Table 2, together with the number of such observations that yielded a positive detection (see also Figs. 1-3), so that the reader can readily grasp the extent of the library. There are several possible outcomes of an observation at the location of a known pulsar: (1) no detection, (2) detection of a point source (for XMM-Newton and Chandra only), (3) detection of pulsation corresponding to the expected pulse period and/or its 11 http://cxc.harvard.edu/toolkit/pimms.jsp 12 https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/ w3nh.pl harmonics. We have enumerated the number of occurrences of each type of outcome in Table 1 for the survey as a whole, and by satellite. There are in total ∼37,000 observations of pulsars, yielding 1599 individual detections in the library, of which 1260 show pulsation and can be expected to yield pulse profiles. The detailed breakdown by individual pulsar is provided in Table 2.
For each known SMC pulsar detection in the XMM-Newton and Chandra observations, we have extracted single source event lists 13 , spectra, source (background and background-subtracted) light-curves, periodograms, and folded light-curves in the broad (0.2-12), soft (0.2-2 keV), and hard (2-12 keV) energy bands. Similar products have been produced for RXTE with a few differences (single energy-band 3-10 keV light curves which are generated on a per-observation rather than per-pulsar basis, no individual source event files). We have combined the three satellite detections together, to explore the long term behavior of each pulsar. Of course, some sources have been detected by one, two or all three satellites.
In the following sections we present an overview of the multi-satellite results: (Section 3.1) presents the distribution of count-rates and L X ; Section 3.2 provides examples of the history and time evolution of properties for individual sources and long-term period derivatives; Section 3.3 shows how pulse period variations can be analyzed for each individual pulsar; and Section 3.4 describes the online catalog and public access to the library. Figure 4 shows the X-ray luminosities from all observations we mentioned and from all three telescopes. The observations cover the range L X = 10 31.2 -10 38 erg s −1 . Filled and unfilled symbols indicate that pulsations were detected or not detected, respectively. We see that RXTE has provided the majority of the detections, since it has observed the most frequently (approximately every week for 15 years). The additional value of Chandra and XMM-Newton is clear in two ways. Firstly these imaging instruments are able to distinguish pulsed emission from un-pulsed. Secondly they are in the form of many detections at lower luminosity, which increases the dynamic range. This is an important attribute because the transition between different accretion modes likely occurs below the typical sensitivity of RXTE. Figure 4 provides valuable information on the disappearance of pulsation. And in this context we have recently performed an analysis of the propellor transition in this ensemble of pulsars (Christodoulou et al. 2016). Minor differences between Christodoulou et al. 2016 and our Fig. 4 are due to inclusion of XMM-MOS data in the former work, and slightly different S/N screening criteria.

Luminosity
In terms of detecting pulsations, a useful empirical result on the sensitivities of XMM-Newton, Chandra and RXTE to pulsars as a function of count rate was obtained by analyzing the library. By constructing histograms of source counts for all positive pulsation detections (s 99%) the distributions of pulsation amplitude, and number of events per detection can both be examined. Figure 5 shows the number of observations with pul-   sations detected by XMM-Newton at s 99%. The abscissa shows the actual EPIC PN photon counts that yielded the pulsation detections in the histogram. The ordinate of the blue histogram is the number of sources from the observations with pulsations detected in each bin (in intervals of 100 counts). The distribution indicates a threshold of about 200 counts below which our ability to detect pulsations becomes rapidly diminished. The distribution of detected pulsations peaks at around 300 counts for XMM-Newton and represents the completeness limit of the survey. Observations should therefore be designed to obtain 200+ counts when searching for periodicities. The decline towards lower net counts is due to the expected decline in ability to resolve pulsations in fainter sources and shorter observations, while the decline toward higher net counts reflects the underlying distribution of pulsar luminosity. The red histogram is the cumulative number of observations with pulsations detected as a function of photon counts. In total, there are ∼70 XMM-Newton observations with pulsations detected in the SMC. Figures 6-8 show similar distributions for pulse amplitude for all three satellites (expressed in net count rate). They are the histograms of all known pulsars detected with s 99% by XMM-Newton, Chandra, and RXTE, respectively. The abscissae are equally binned in logarithmic count rate. It is apparent that a turnover occurs just above log 10 A −1.2 or 0.06 count/s for both ACIS and PN, and for the RXTE amplitude histogram a simi-lar turnover occurs at log 10 A −1.7 or 0.02 count/s. To compute the completeness fraction of the XMM-Newton, Chandra and RXTE surveys the underlying pulsar L X and P F distributions would need to be obtained first. We hope that modeling efforts will be motivated by statistical results such as these.

Time Evolution of Source Properties
In the combined pipeline, the standard data products for each pulsar include the pulse period, luminosity, pulsed fraction, amplitude, and the spin period detection significance. The time evolution of these quantities yields the period derivativeṖ , the accretion torque, the orbital period, and the duty cycle of each source. Time-series plots of these quantities are provided as multi-panel PDF figures in the online journal for all pulsars in the sample, and the underlying data are provided in machine readable form. These time series can be read as a history of the on/off status and outburst state of each pulsar as a function of time. Such information is useful for comparison with other time-domain databases extending over the same time period, for example the OGLE and MACHO optical monitoring facilities. The plots are also intended as a useful resource for researchers interested in selecting their own sub-samples for further analysis or modeling.
As an example, the results for the HMXB pulsars SXP348 and SXP1323 are illustrated in Figs. 9 and 10, respectively. Each plot has 5 panels, with the data from each satellite plotted with a different color and symbol.
The first (top) panel is the time series of L X in units of erg s −1 following the cross calibration scheme described in section 2.4.2. These are in fact the same points that appear in Fig. 4. L X values are plotted for all positive point-source detections (whether pulsed or not) from Chandra and XMM-Newton, and for RXTE pulsed detections.
The second panel reports the pulsed fraction for Chandra and XMM-Newton detections with pulsations at s 99% with solid symbols. Detections at s < 99% are open symbols.
The third panel shows pulse amplitude in units of count rate, where the Chandra and XMM-Newton rates are in units of PN count s −1 after scaling as described in Section 2.4.2, and RXTE rate is in units of count PCU −1 s −1 . We have followed the practice of earlier works (e.g. Galache et al. 2008) by plotting RXTE values even in cases of no pulsation detection; these open points are intended to represent upper limits.
The fourth panel shows the spin period (P ), with values plotted only for cases of detection at s 99%. A linear fit is also displayed in order to highlight the long-term trend in period derivative. In the example of SXP348 (Fig. 9) two epochs of spin-up are seen, each lasting about 500 days, between which the pulsar returns to its long-term average value. In our second example, SXP1323 (Fig. 10) there is an overall long-term spin-up trend, with other apparently organized variations superimposed. This analysis was performed for all the pulsars with sufficient data to do so, and the slopes of these linear fits, their uncertainty, and the standard deviation of the points around the fit are collected together in Table 3. Here standard deviation implies how much the period of the pulsar varies on long timescales: a low standard deviation indicates that the data points tend to be close to the best fitting line of the set, while a high standard deviation indicates that the data points are spread out over a wider range of values around the best fitting line.
The fifth (bottom) panel reports the significance of the highest power-spectrum peak in the search-range for that pulsar (See section 2.4.1). Pulsation searches were only performed for light curves with greater than 50 counts in the case of XMM-Newton and Chandra. The solid symbols denote s 99% and in turn dictate how the points in other panels are plotted.

Long-term Period Variations
As expected for accreting pulsars in binary systems, all the pulsars in our library show variations in their pulse periods. These variations stem from periodic doppler shifting of the frequency due to orbital motion and from the exerted accretion torques. The first type of peri-odic variation may reveal the orbital period of the binary (see for example Townsend et al. 2013), whereas the second type of secular variation gives us a window into the physics of angular momentum transfer and dynamical evolution in binary systems (Davidson & Ostriker 1973;Ghosh & Lamb 1979). The two types of variations can be seen in some of the multi-panel figures included in the online journal as supplementary material. The orbital modulation in some sources is not very clear, e.g., SXP 348 in Fig. 9, due to the few spin period data points. Schmidtke & Cowley (2006) found weak evidence for an orbital period of 93.9 days for SXP 348 in OGLE-II data. However, Schmidtke et al. (2013) reported from re-analysis of the same data (plus OGLE-III and -IV data in which it was not found), that the period was an artifact. Coe & Kirk (2015) simply quoted its orbital period as 93.9 days. From MJD 51000 until MJD 52000 is about 10 orbital cycles (with an orbital period of 94 days). During MJD 51000-52000, it only shows spin down and up once. SXP1323 (Fig. 10) has no secure orbital period known according to Coe & Kirk (2015).
Here we carry out a linear regression of the spin periods of 53 pulsars for which we have enough data in order to search for secular drifts on timescales much longer than the orbital periods.
We have investigated the entire spin-period history of each pulsar in our library by calculating the best-fit slope to the measured spin periods and its corresponding errors, as described in Section 3.2. The best-fit slopes, i.e., the measuredṖ values, are listed in Table 3 and they are also included in the multi-panel history plots such as Figs. 9 and 10. At the = 1.5 level or better, we find that 12 pulsars spin up and 11 pulsars spin down, while 30 pulsars appear to haveṖ ≈ 0. The "C" in Table 3 does not mean that the 30 pulsars do not go through spin-up and/or spin-down episodes, but that no net changes of significance are observed over the ∼15 year duration of the survey. Some sources could well have undergone spin period changes that have averaged out over the survey period. In our dataset, 5 pulsars only have one single detection with s 99% pulsations found as shown in Table 3. The remaining 9 pulsars do not have significant pulsations detected at all.
The cumulative outcome of this analysis is illustrated in Fig. 11 in which we plot the known orbital periods P orb versus the known spin periods P (Coe & Kirk 2015) with additional color-coded information about the sign ofṖ . Green and red symbols indicate that the pulsars spin up (Ṗ < 0) or down (Ṗ > 0), respectively (e.g., SXP1323, Fig. 10). On the other hand, blue symbols indicate pulsars have noṖ value due to lack of observations or lack of pulsations detected in our survery. Unfilled symbols at the bottom of Fig. 11 show the pulsars for which P orb is unknown; currently, 43 SMC pulsars in our library do have measured orbital periods.
The relation between the orbital and spin periods of 43 pulsars in Fig. 11 can be reasonably matched by a power law. The best-fit power law to the data is described by the equation  This equation shows that the longer the orbital period the wider the binary orbit, hence lower mass transfer with lower specific angular momentum occurs, and ultimately the pulsar rotation ends up being slower. Fig. 11 is effectively our updated version of the well-known Corbet (1984) diagram for the SMC. Overall, the HMXB spin period distribution has been shown to be bimodal by Knigge et al. (2011). Those authors suggested that there are two underlying populations tesulting from the two distinct types supernovae-electron capture and ironcore-collpase-that produce neutron stars.
3.4. Public Access to the Library We will incrementally release this library to the public, from which one can download all of the light curves, periodograms, point source event lists, spectra (within the soft, hard and broad energy bands) of each known pulsar in the SMC & LMC. For this paper, release 1.0 provides the long-term observable parameters (i.e., count rate, luminosity, spin period, pulsation amplitude, pulsed fraction) of the known pulsars in the SMC with the combination of all 3 satellite detections. The library is provided in catalog form and as a series of multi-panel PDF figures (58 in total).
Note. -Period derivativeṖ and its error are determined from a linear fit to the time series of P spin measurements. σs is the standard deviation indicating how much the period of the pulsar varies on long timescales around the linear fit. is the ratio ofṖ to its own error as described in Section 3.3. Classification is based on the sign and value of : U: −1.5 (long-term spin up), D: +1.5 (long-term spin down), C: −1.5 < < +1.5 (consistent with no long-term variation). In the last column, # of points is the number of detected pulsations used in the linear fit.   Table 3. Blue symbols represent pulsars whose spin period derivatives are unknown due to lack of observations or lack of pulsations detected.
are generally needed for reliable pulse detection, as well as subsequent spectral energy analysis at moderate resolutions.
The library spans ∼15 years and includes all key observational parameters generated by individual pipelines for the 3 X-ray telescopes which, in turn, are combined to produce a multivariate time-series for each X-ray pulsar (as illustrated in Figs. 9 and 10). These time series are provided in PDF and/or machine readable table form for each of the 67 SMC pulsars.
One of these physical parameters, luminosity, as the function of the spin periods can map out the boundary of the pulsars in the Small Magellanic Cloud as 10 31.2 to 10 38 erg s −1 . The fast pulsars are less frequently observed. They are rarely experiencing outburst. Once they are in the outburst state, they are more likely to be the giant outburst. The very few detected fast pulsars tend to have a higher luminosities with XMM-Newton observations. They might go through the outburst state, e.g., SMC X-2 (Palombara et al. 2016), SMC X-3 (Townsend et al. 2017;Weng et al. 2017). From the histogram Figs. 1-3, by chance, the long spin period pulsars are more frequently observed. The imbalance may guide us to propose the observing proposals for the fast pulsars.
For this paper introducing the library, we have performed an analysis of long-term period derivatives as an example of the types of investigation that can be done with this large sample of pulsars. The data show that long-lived accretion torques are present in about half of the sample. Okazaki et al. (2013) have pointed out that the standard picture of spin-up occurring only during transient outbursts is difficult to reconcile with the fact that the viscous timescale in the accretion disks is often longer than the orbital period. Nonetheless, the transient spin-up paradigm grew out of observations of the few pulsars with long-term period monitoring such as EXO 2030+375 (Wilson et al. 2008). By fitting the long-term pulse period values with linear models, we have determined that, at a level of = 1.5 or better, 12 SMC pulsars are spinning up and 11 spinning down on the ∼15 year timescale of the survey. On the other hand, 30 pulsars are consistent with no significant long-term average spin period changes and the remaining 14 pulsars remain uncategorized due to paucity of high signalto-noise spin period measurements. The largest longterm average spin-up and spin-down values are -6.5e-03 ± 2.8e-03 s/day (SXP1323) and 7.1e-03 ± 2.2e-03 s/day (SXP1062), respectively. For comparison, such a very long-lived spin-up has been previously observed in just a few systems (e.g., Krivonos et al. 2015). The relationship between the long-term spin-up rate and the X-ray luminosity has been the subject of a recent study  using the Fermi γ-ray observatory and the MAXI X-ray monitor on the International Space Station (ISS). This study found thatṖ and L X are closely correlated (−Ṗ ∝ L X 6/7 ), as was also found by Coe et al. (2010) for a sample of SMC pulsars. With this larger set ofṖ values it will be possible to explore variations in the relationship predicted by Ghosh & Lamb (1979).
The initial public release of our library that accompanies this paper includes the long-term observable properties of the known pulsars in the SMC. Future releases will include mid-and high-level data products for every individual pulsar detection. These products will include high time-resolution (raw and folded) light curves, periodograms, single source event lists, and calibrated energy spectra. We are also working to add to the library the 14 Be/X-ray pulsars in the LMC. As new observations enter the archive, our library will be updated and it will serve as a permanent resource for the community. Furthermore, we note other contemporary data-mining efforts that are under way to discover pulsars in the archival data (Israel et al. 2016).
The ultimate goal of this project is to unleash the power of statistics on a large, unbiassed, observational sample of Magellanic Be/HMXB pulsars. These pulsars constitute precisely such a sample due to the presence of a large population of HMXBs within a small volume of space, at a known distance, that are also embedded in a readily resolvable stellar population. The rich statistical results from our library will certainly improve our understanding of magnetized accretion processes in HMXBs as the various products (pulse profiles, periodograms, X-ray spectra) will be used in the testing of theoretical models of gas flows through NS magnetospheres.