A Search for Pulsars in Ultracompact X-ray Binary Systems

MASTER'S THESIS

A Search for Pulsars in Ultracompact X-ray Binary Systems

Bastian Hacker 2013

Master of Science (120 credits) Space Engineering - Space Master

Luleå University of Technology

Department of Computer Science, Electrical and Space Engineering

A search for pulsars in

ultracompact X-ray binary systems

Master Thesis in astrophysics

Bastian Hacker 06/06/2013

supervisor: Matteo Bachetti advisor: Natalie Webb

Galaxies, High Energy Astrophysics and Cosmology (GAHEC) group, Institut de Recherche en Astrophysique et Planétologie,

Luleå tekniska universitet,

Université Paul Sabatier Toulouse III

Abstract

We present a search for periodic pulsations in the X-ray emission of seven Ultracompact X-ray Binary (UCXB) systems, in which the coherent time domain resampling method was used to recover signals distorted by the target’s orbital motion. While other UCXBs have been observed as X-ray pulsars, our targets still lack such a detection. We demonstrated the superior sensitivity of our search method over standard approaches, enabling us to find very faint signals at the cost of a large computing effort. However no pulsations were detected. Therefore we deduced stringent upper limits for the pulsed fraction in each observation, which are on the order of 0.5%.

Acknowledgements

This work was carried out at IRAP, the astrophysics and planetology research institute in Toulouse, France as a final thesis in both the international SpaceMaster programme and the ASEP astrophysics master course. The possibility to conduct the project in this professional and inspiring environment is kindly appreciated.

I want to thank my supervisor Matteo Bachetti for his invaluable support throughout the project. It was a real pleasure to work together and learn so much about astrophysics and computing over countless productive coffee sessions. Thank you Matteo for having been an excellent supervisor!

Thanks also to my advisor Natalie Webb for her professional support in guiding the project, her precious advice and resources. Her assistance helped to improve the work substantially and was very appreciated.

I’d like to thank Luciano Burderi, Tiziana di Salvo and Alessandro Riggio who laid the groundwork to apply the coherent pulsar search technique together with my supervisor Matteo Bachetti. It was their ideas that created the project and their previous work and experience that I could build on.

I also want to express my gratitude to all the people who organised and carried out the study course, especially Soucail Geneviève who assisted me throughout my academic period in Toulouse and gave precious guidance.

I’d like to thank my research group members for the fruitful discussions and my student colleagues for their valuable advice and distraction when needed.

Last but not least I thank the CINECA supercomputing centre for granting nearly two million computing hours on their Fermi supercomputer. The project would not have been possible otherwise.

Contents

1 Introduction 4

1.1 Pulsars . . . 4

1.2 X-ray emission from pulsars and its detection . . . 5

1.3 X-ray binaries and accreting millisecond pulsars . . . 6

1.4 The pulse signature of a binary pulsar . . . 8

1.5 Search techniques . . . 8

2 Methods 10 2.1 The xpsrspec software . . . 11

2.2 Supercomputing . . . 12

2.3 Time domain resampling . . . 12

2.4 Residual timing errors . . . 13

2.4.1 Eccentricity . . . 13

2.4.2 Post-Keplerian Parameters . . . 13

2.4.3 Varying pulsation period . . . 14

2.4.4 Positioning error . . . 14

2.5 Search grid . . . 15

2.6 Parameter constraints . . . 18

2.7 Fourier analysis . . . 20

2.8 Follow-up epoch folding . . . 22

3 Analysis 22 3.1 Tests with a known pulsar . . . 22

3.2 Candidate binaries . . . 24

3.3 Results . . . 25

4 Conclusions 27

References 28

1 Introduction

This work presents a search for periodic pulsations in the X-ray emission of Ultracompact X-ray Binary (UCXB) systems. These systems may contain a pulsar whose signal gets smeared out due to the binary orbital motion, which makes the detection difficult. We implemented a program that recovers possible pulsation signals and used it to search seven of the most promising UCXB systems.

In this section I introduce all the key concepts, in section 2 I describe the methods and steps I used in the search and in section 3 I report a test of our tools on one system with known pulsations and the results of our search.

1.1 Pulsars

Pulsars (pulsating stars) are neutron stars (NSs) emitting fast pulsations. They were first discovered serendipitously in 1967 through radio observations (Hewish et al., 1968). Soon after their discovery it was established that pulsars must be rapidly rotating, highly magnetised NSs (Gold, 1968; Pacini, 1968). They have theoretical masses of the order of 1.4M (solar masses) packed into a very small sphere with a radius of around 10 km (Lattimer and Prakash, 2007).

Their density is higher than inside the atomic nuclei. This makes NSs some of the most extreme objects in the universe and a unique celestial laboratory for fundamental physics.

Increasing numbers of pulsars (now about 2000) have been detected, most of them within our galaxy. The majority spin at periods of the order of a second while about 10%, the so called millisecond pulsars (MSPs), spin with periods down to few milliseconds (figure 1). Much can be learned from the distribution of spin periods and their derivative in time, but for this work it is sufficient to notice that MSPs are more likely found in a binary system (about two thirds of them).

10

10

10

10

10

P [s]

10

10

10

10

10

10

10

10

10

10

10

10

dP/dt [s/s]

radio pulsars

high energy pulsation binary pulsars

Figure 1: Distribution of pulsation periods P and spin-down rate of 1902 known pulsars from the ATNF catalogue (Manchester et al., 2005). The two distinct groups of “normal” and millisecond pulsars can be identified along with a tendency of millisecond pulsars to be found in a binary system.

Pulsars can have extremely stable pulsation signals that are in some cases comparable to atomic clocks. This is used to study the pulsar rotation behaviour precisely, determine the pulsar position and movement very accurately or even for precision tests of general relativity (Taylor and Weisberg, 1982; Lyne et al., 2004). Therefore pulsars are precious tools in astrophysics and the search for new and different exemplars is a worthwhile endeavour.

1.2 X-ray emission from pulsars and its detection

Pulsars carry very strong magnetic fields (illustrated in figure 2, right) and generally lose rotational energy by very low-frequency magnetic dipole emission into their surrounding medium which is not observable directly. The pulsed emission of radio waves is produced by charged particles that are accelerated along the magnetic field lines and therefore emit curvature and synchrotron radiation. X-ray pulsed emission on the contrary happens when matter is accreted onto the pulsar’s surface. The infalling particles lose so much gravitational energy while they are magnetically channeled onto the magnetic poles that they form a hotspot of several million degrees on the surface which emits strong thermal X-rays (White et al., 1983). A directed beam is then formed by the opaque accretion flow, which absorbs most escaping photons except the ones near the magnetic pole direction. Both mechanisms create an apparent pulsation with the frequency of the pulsar’s rotation when the magnetic poles don’t conincide with the rotational axis and the rotating light beam grazes the line of sight towards the observer.

Most pulsars are seen and detected by pulsed radio emission with large-area radio telescopes.

Only a fraction of these is also visible in high energy bands such as X-rays (marked in figure 1 with crosses). On the contrary X-ray pulsars in X-ray binary systems – which are being targeted in this work – have so far been invisible in the radio domain (Lorimer, 2008) with only one recently discovered exception (Papitto et al., 2013). The dense ionized clouds around these pulsars screen the radio radiation effectively while X-rays can escape. Therefore X-ray observations are necessary to detect these pulsars and learn about them during their close binary phase.

This work deals purely with X-ray signals. The X-ray radiation is detected as discrete photons at individual arrival times. All the available data is collected by satellites because the earth’s atmosphere is opaque to X-ray radiation. But even above the atmosphere one requires large collecting areas and long exposures to reach reasonable signal levels. Additionally the search for high frequency pulsations requires very accurate timing. Both of these are provided by NASA’s

Figure 2: Left: Artist’s illustration of a LMXB system with the primary in blue, the secondary in yellow and an accretion disc in red. Right: Illustration of an accreting pulsar. (Credit: NASA)

Rossi X-ray Timing Explorer (RXTE) satellite¹, which has been in operation from 1996 to 2012.

Its collected observation data on the HEASARC archive² has been the basis for numerous pulsar searches and it is also the data source for this work.

RXTE was operated in a low-earth circular orbit at 580 km altitude with a 90 min orbital period, where observations can get interrupted by earth occultations. Of its several detectors we use the PCA (Proportional Counter Array) because of its large collecting area, good time resolution and directional selectivity. It is sensitive to X-rays of energies between 2–60 keV with 6500 cm² collecting area and has a spatial resolution of 1^◦ FWHM thanks to a collimator.

This relatively coarse spatial resolution comes as an issue for crowded fields such as globular clusters which contain many of the otherwise possible candidates for this project, because the signal of each source therein is superimposed with numerous other emissions that act as noise background. The time resolution is 1 µs for observations in the “Good Xenon” mode and 122 µs in the “Science Event” mode, which are both well-suited for accurate timing measurements even for the fastest expected millisecond pulsations. The dead time of the PCA instrument is varying for different possible settings between 12 µs and 50 µs. This is around two orders of magnitudes less than the typical count rates of 100 evt/s so its effect on the statistics will be negligible.

1.3 X-ray binaries and accreting millisecond pulsars

Binary stars are gravitationally bound systems of two stars orbiting each other. They constitute a significant fraction of all star systems in the galaxy. Their orbital motion around the common center of mass obeys Kepler’s laws, in particular

4π²

P_orb² = G(M_p+ M_c)

(a_p+ a_c)³ (1)

where G is the gravitational constant, P_orb the common orbital period, a_p and a_cthe semimajor axes and M_p and M_c the masses of the primary and the companion star respectively. When the system is viewed at an inclination i (0^◦ being face-on) one can often only measure the projected semimajor axis of the primary along the line of sight X = a_psin i and P_orb, from which one can infer the so-called mass function f_X by modifying eq. (1) with the center of mass equation Mpac= M_cap:

fX ≡ 4π²X³

G P_orb² = (M_csin i)³

(M_p+ M_c)² (2)

The mass function gives a combined constraint on the masses M_p and M_c which is determined by measuring X and P_orb, while additional knowledge is required to infer each mass individually.

It will be useful when we are going to determine parameter constraints for our investigated systems.

Binary stars are subject to stellar evolution over time much like single stars, ending their main sequence lifetime sooner the more massive they are. Our systems of interest start out with a primary between about 10M and 25M (Heger et al., 2003), which will cease in a supernova explosion before its companion, leaving behind a neutron star. In case the binary is not disrupted in the explosion, the primary will subsequently start accreting mass from the companion via stellar wind from the companion (Bondi-Hoyle accretion) or Roche lobe overflow (when the companion surface expands beyond the inner Lagrangian point, enabling matter to flow towards the primary, illustrated in figure 2, left), gaining so much energy from the infalling matter that the binary emits bright X-ray radiation: an X-ray binary.

1http://heasarc.gsfc.nasa.gov/docs/xte/XTE.html

2http://heasarc.gsfc.nasa.gov/docs/archive.html

Table 1: The 14 first known AMXPs with orbital parameters and root mean square pulsed fraction (pf)

Name Found fpulse[Hz] Porb[s] X [ltms] pf Ref.

SAX J1808.4–3658 1998 400.975 7249.16 62.81 4.1%–7% ¹

XTE J1751–305 2002 435.318 2545.34 10.11 3.5% ²

XTE J0929–314 2002 185.105 2614.75 6.29 3%–7% ³

XTE J807–294 2003 190.624 2404.42 4.82 3.1%–6% ⁴

XTE J1814–338 2003 314.356 15388.72 390.63 12% ⁵

IGR J00291+5934 2005 598.892 8844.09 64.99 8% ⁶

HETE J1900.1-2455 2007 377.296 4995.26 18.41 1%–3% ⁷ SWIFT J1756.9-2508 2007 182.066 3282.10 5.94 4.2% ⁸

Aql X–1 2008 550.27 68213 N/A 1%–10% ⁹

SAX J1748.9-2021 2008 442.36 31550 390 1%–3% ¹⁰

IGR J17511-3057 2009 244.834 12487.51 275.20 6%–9% ¹¹

NGC6440 X-2 2009 205.892 3438 6.22 7.5% ¹²

Swift J1749.4-2807 2011 517.920 31740.73 1899.53 6%–23% ¹³

J17498-2921 2011 400.990 13835.62 365.11 10%–20% ¹⁴

1Wijnands and van der Klis (1998),²Markwardt et al. (2002),³Galloway and Chakrabarty (2002),⁴Markwardt et al. (2003),⁵Markwardt and Swank (2003),⁶Galloway et al. (2005),⁷Galloway et al. (2007),⁸Krimm and Markwardt (2007),⁹Casella et al. (2008),¹⁰Altamirano et al. (2008a),¹¹Markwardt et al. (2009),¹²Altamirano

et al. (2010),¹³Altamirano et al. (2011),¹⁴Markwardt (2011)

X-ray binaries are classified in two categories: High-mass X-ray binaries (HMXBs) and low-mass X-ray binaries (LMXBs) depending on the companion (see Verbunt, 1993, for a comprehensive overview). In HMXBs the companion has a mass greater than 1M and usually consists of an early-type main sequence star which makes the binary appear optically bright.

HMXB systems are not at the center of interest of this work, because they are supposedly young which makes them unlikely to harbour millisecond pulsars that we are searching for.

In LMXBs on the contrary the compact primary has a < 1M companion, which is either a late-type main sequence star or a white dwarf. The primary in LMXBs is always a neutron star or a black hole. Since the companion is also relatively small, LMXBs are optically faint and the X-ray emission from the accretion will dominate the flux. The companion is also too small to evolve further than the white dwarf state, which allows for a long term (10⁹yr) mass transfer onto the primary to speed up its rotation with angular momentum from the binary motion.

UCXBs (see Nelemans and Jonker, 2010, for a recent review) are the most extreme case of LMXBs with an orbital period below about 80 min, which by eq. (1) also requires very tight orbits (. light seconds, lts). UCXBs are particularly interesting as their tight orbits cause measurable angular momentum losses by gravitational wave emission (Hulse and Taylor, 1975).

This causes the binary distance to shrink together with the Roche lobe of the companion, which generates continuing mass and angular momentum transfer onto the primary.

The long-term accretion onto a neutron star primary can spin it up to very high frequencies by transferring orbital angular momentum onto the NS (as for instance shown for IGR J17511-3057, Riggio et al., 2011), potentially leading to a millisecond pulsar. According to the recycling scenario (Bhattacharya and van den Heuvel, 1991) former binary pulsars, which switched off their pulsation after their magnetic field decayed, may be activated again during the accretion process and emerge as a recycled millisecond radio pulsar after the accretion has finished (Archibald et al., 2009). During the accretion process the radio pulsation is screened by the accreting matter (Cumming et al., 2001), so that accreting millisecond radio pulsars can hardly be detected.

However pulsed X-ray emission can be visible. To date 14 of such Accreting Millisecond X-ray Pulsars (AMXPs) have been found (Patruno and Watts, 2012), starting with the discovery of pulsations from SAX J1808.4-3658 (Wijnands and van der Klis, 1998). We list them in table 1 to have a first estimate on the distribution of parameters in such systems. All AMXPs are transient systems with periods of quiescence (low luminosity) and outbursts (high luminosity) which correspond to a varying accretion rate onto the neutron star. Furthermore some systems are intermittent, which means that pulsations were not always detectable even during outbursts.

The question which factors exactly cause an accreting NS to exhibit pulsations or not is still a matter of active discussion (Özel, 2009) and may be advanced by further detections.

1.4 The pulse signature of a binary pulsar

C

observer P

Figure 3: Sketch of a LMXB system with a NS primary (P) and low mass Roche lobe filling companion (C) orbiting each other. Pulsations from the NS are received modulated by the orbital motion.

UCXB sources are often far away (∼ kpc, see e.g. review in Lorimer, 2008), appear faint (. 1 photon/pulseperiod) and have only a low fraction of their emission pulsed (pulsed fraction, pf, see table 1). Therefore single pulsations cannot be seen in the data, but only averaged pulse profiles after thousands of pulse periods. This averaging only leads to the right profile if the behaviour doesn’t change over time or time variations can be corrected for.

What originates from a binary pulsar as strictly periodic pulsation in its co-moving reference frame at time t_em may not be perceived as such by an outside observer at t_rec. The orbital movement of the pulsar along the line of sight towards the observer d(t_em) introduces a varying time delay due to the finite propagation speed of light c:

t_rec = t_em+ c · d(t_em) (3)

This goes along with a varying (nonrelativistic, v . 10⁻⁴c) Doppler shift of the pulsation frequency. The pulsation signal received by the observer is therefore frequency-modulated (figure 3) and will show up with a significantly lower maximum amplitude in a simple Fourier power spectrum. This distortion prevents a straightforward detection of the signal in many cases, especially in the most extreme and therefore interesting systems. A variety of countermeasures has been developed consequently, one of them being the subject of this work.

1.5 Search techniques

Finding periodic modulations in any sort of signal is one of the most classic tasks in science and engineering. The one most important technique for this is the Fourier Transform, because in

its implementation as discrete Fast Fourier Transform (FFT) it can analyse all frequencies in a given evenly spaced signal of length n very efficiently with a time complexity of O(n log(n)).

Therefore it has become one of the most standard techniques in X-ray timing for unknown pulsed emissions (Leahy et al., 1983; Van der Klis, 1988).

Other techniques apart from the FFT have also been applied in the search for periodic emission, namely the Lomb-Scargle periodogram (Wen et al., 2006) and epoch folding. These bring advantages in sensitivity in some cases, are however more computationally intensive and don’t alleviate the problem of frequency modulation.

The recovery of binary modulated signals in particular has been investigated for more than two decades. Viewed on observation time scales shorter than a small fraction of the binary orbital period, the Doppler shift is almost constant and a single delta peak in the signal power spectrum can be observed, provided the source emits enough signal photons during this time frame (Wijnands and van der Klis, 1998). This task is performed easily and efficiently on large data archives and therefore we have to assume that all strong and clear pulsation signals in our data are already discovered by now. Similarly if a long observation signal is Fourier transformed without considering the modulation, strong pulsations remain detectable as a pattern around the pulsation frequency with reduced amplitude (Jouteux et al., 2002).

An approach that is often taken is the dynamical power spectrum (Lorimer, Duncan Ross and Kramer, 2005), which stacks many power spectra of subsequent signal slices in a time–frequency diagram. Frequency-drifting signals which are barely significant in individual subspectra may then be identified as traces in the 2D diagram. Several methods to exploit this were developed such as the application of the Hough transform (Aulbert, 2006).

A sensitivity improved technique is the acceleration search (also called quadratic coherence recovery), where the signal is corrected for several constant trial frequency drift rates prior to the Fourier Transform. While this allows to fully compensate for small intrinsic frequency drifts on astronomical timescales, application in binary systems (Wood et al., 1991) still only allows to recover observation lengths of a fraction of the binary period and can yield a limited improvement in sensitivity over the uncorrected search on the order of 2.5 (Johnston and Kulkarni, 1991).

This technique however cannot exploit available observations which are much longer than the shortest orbital periods.

The straightforward approach to fully recover the pulsar signal is to correct the photon arrival times (time domain resampling) for the exact binary motion and then Fourier transform the result. This is also called orbital-parameter-fitting coherent search because it recovers the signal’s coherency over whole observations. In this approach the space of a priori unknown Keplerian orbital parameter sets is searched by brute-force. Unfortunately the necessary number of parameter sets for this is typically large and made the search infeasible in the past, although its potential power was recognized long ago (Van der Klis, 1988). But ever growing computing power alleviates the problem so that groups like Dib et al. (2005) managed to search coherently for pulsations in the LMXB 4U 1820-30, where P_orb is known very precisely. Further work like by Bachetti (2010) demonstrated that coherent searches can be performed with even coarser parameter constraints by now. Both these searches found no new pulsations despite their superior sensitivity. In addition coherent methods are routinely applied in computationally intensive binary radio pulsar searches for certain ranges of orbital periods (Allen et al., 2013). The work of Bachetti (2010) constitutes the footing on which this project builds, providing software tools and a theoretical framework which is used here and extended at some points.

At the same time approaches to mitigate the vast computing effort were made. Most notably the sideband (or phase modulation) search technique developed by Jouteux et al. (2002) and Ransom et al. (2003). They exploited the fact that the theoretical power spectrum of an

uncorrected binary pulsar signal observed over several orbital periods is not simply a smeared- out peak, but has a regular comb-like pattern due to the signal’s very coherent and periodic nature. In essence they proposed to take a second FFT on subsequent parts of the raw signal’s power spectrum to convert the hidden regular pattern into one large detectable peak. This is computationally efficient as it requires no free parameter adjustment. The downside of the sideband search is that it reduces the detectable signal power by a factor of 0.4 compared to the coherent search even in the best case simulations and drops further if the number of orbits per observation is small or X · f_pulse is high. This places the method between the primitive uncorrected searches and expensive coherent searches.

Another noteworthy way to speed up the calculation is the time differencing method (Atwood et al., 2006) originating from the gamma ray domain. Since gamma ray count rates are low and therefore long observations are required, the huge resulting size of the FFT array is a particular problem there. Their cure is to collect all differences of corrected photon arrival times up to a certain maximum value and take the power spectrum from that instead of the complete time series. While the method performs well for their application it does no good for the X-ray analysis of LMXB pulsars. Due to the much higher count rate of X-rays the effort savings from shorter FFT arrays are outweighted by the summing process if the maximum time difference allowed is larger than on the order of a second. But the sensitivity of the method decreases with smaller time differences such that in our case at least an order of magnitude would be lost, which would undo all the advantages of parameter correction.

Therefore the classic coherent time domain resampling search yields the highest sensitivity for close binary systems that we investigate. The absence of such reported high-sensitivity searches for many of the known candidate systems leaves an opportunity to use available computing power and systematically work through the latest observation archives. The results can either be the discovery of previously hidden weak signals that open the way for further investigation or tighter limits on possible pulsation amplitudes, which supports the discussion of why pulsations are suppressed in some of the systems.

2 Methods

In this work we take RXTE observations of candidate UCXB systems from the Ritter catalogue (Ritter and Kolb, 2003, update RKcat7.19, 2013) to search for previously undiscovered X-

ray pulsations. Our method of choice is the computationally-intensive coherent time domain resampling method (sec. 2.3) with FFT analysis for peaks in the power spectrum (sec. 2.7).

After this first step to retrieve candidate parameter sets we run a follow-up epoch folding test (sec. 2.8) on the candidates to verify or reject them with greater fidelity.

Pulsars in binary systems may be overlooked with fast standard search techniques and remain buried in observation archives until analysed with greater sensitivity (e.g. like in Mickaliger et al., 2012). Our search method reaches the greatest sensitivity of all search methods in sec. 1.5 and hence there are chances to find weak signals which have been overlooked by previous less sensitive searches. My search was performed as follows:

Selecting appropriate candidate UCXB systems from catalogues, which have an estimate of the orbital period but no known coherent pulsations. We are only interested in the most compact systems because these are the ones for which the coherent search can provide a better sensitivity over other more computational efficient techniques (see sec. 1.5).

Collecting previously known constraints of the candidates from the literature such as the orbital period, orbit inclination, companion type and previously observed burst oscillations

(brightness variations from nuclear burning material rotating close to the NS’s surface, Strohmayer et al. 1996) and determined all parameter constraints based on the discussion in sec. 2.6.

Data reduction: I downloaded RXTE observation data sets from NASA’s public HEASARC archive and chose suited ones with large observation duty cycles. I preferred observations with T < 3.5 hr because longer ones would not fit in the computer memory at a time (see sec. 2.2). I refrained to select observations by spectral states which may be connected to the occurrence of pulsations (Altamirano et al., 2008b), avoiding the danger of premature rejection. The data are available as preprocessed event lists, so the only reduction I had to do was to barycenter the arrival times and concatenate separate files of one observation.

Execution of the orbital parameter fitting FFT search (sec. 2.7) on the supercomputer and collected candidate pulsations. For a possible detection the maximum Fourier power has to exceed the theoretical noise level.

Evaluation of the best detection candidates with epoch folding (sec. 2.8) to generate pulse profiles. I used statistical tests to find the profiles most unlikely generated from noise and use them to constrain the maximum possible pulsation amplitude present in the data.

2.1 The xpsrspec software

The main tool for the analysis is xpsrspec, a custom program written in C by my supervisor Matteo Bachetti (Bachetti, 2010) specifically for the task of coherently searching for X-ray pulsations. In his work he used it to analyse a 1 hr observation data set of the UCXB candidate system and X-ray burst source M15-X1 and found that no pulsations above a pulsed fraction of 1% were present. It should be noted that xpsrspec was not an out of the box program that would automatically process a number of given observations; In fact large parts of my work were devoted to make the code work reliably in the supercomputing environment with new and larger data sets. During this project I enhanced xpsrspec’s execution speed by improving the file access scheme and implementing RAM memory savings. I also added additional features such as plotting parts of the power spectrum for diagnosis and using the H-test (see sec. 2.8) to identify the best pulse profiles.

xpsrspec opens observation data event lists in the FITS format and creates a grid of trial parameter sets (sec. 2.5) in the given range (sec. 2.6). The trial parameter sets are distributed to many compute nodes via Message Passing Interface (MPI). Each compute node will then work through its given parameter sets performing the following tasks:

• Time domain resampling of the observation data according to the given parameter set and binning into a regular time grid

• FFT of the time array using the fast fftw library to get a discrete power spectrum

• Searching for peaks in the power spectrum that exceed the noise threshold and saving all associated parameters as a candidate pulsation

After finishing all data sets xpsrspec has an epoch folding mode to extract significant pulse profiles from the collected results and estimate pulsed fraction limits.

We needed high computing efficiency as the total computing effort is large and supercomputing time is precious. The bulk computing effort for typical data sets should be spent for the FFT, as it is the only operation that scales super-linear with the data size. I verified this with the profiling tool Valgrind which showed that over 50% of the time was spent with the FFT for a test data set, so there is no significant margin left for optimisation. The other potentially

critical operation is the write process of candidates to the disk, because large amounts of data from parallel processes must be written at the same time. This task was accomplished using the the Network Common Data Form (NetCDF) library which is designed for high performance binary data access. The write access of xpsrspec is designed to write many candidates at once, so the write frequency could be kept low without slowing down the computation.

2.2 Supercomputing

The large number of trial parameter sets (of the order of 10⁶ for each observation data set) that is needed to perform the coherent time domain resampling search requires huge computing resources, so that the processing on one cpu would take hundreds of years. This may explain why the available data on some candidate systems have not been analysed coherently before. But as all parameter sets can be processed independently the problem is “embarrassingly parallel”, making it well suited for a supercomputer.

The real necessity of high power computing for this problem is for instance illustrated by the Einstein@Home search (Allen et al., 2013), which searched radio, gravitational wave and gamma ray data. Compared to our search they even had an additional degree of freedom, the dispersion measure, which is an a priori unknown radio-frequency dependent time delay of the signal which is introduced by free electrons in the interstellar medium. Consequently their computing effort was even greater and more than 10⁹ cpu hours have been provided by a large number of individuals on their private computers.

In our case we applied for computing resources on the Fermi supercomputer at the CINECA supercomputing centre in Italy which is currently ranked the 9th most powerful machine on the Top500 list³. Our proposal for a “Class C2 project” was accepted in March. This is considered a small project with scientific characterisation and provides 1 966 080 single-cpu hours.

Fermi is a IBM Blue Gene/Q massively parallel computer with each cpu clocked at 1.6 GHz.

The system allows for assignments of at least 1024 parallel computing cores for each program execution. Sub-units of 16 computing cores share a 16 GB RAM memory, which puts a tight 1 GB memory constraint on each thread in order to make use of all available cores.

2.3 Time domain resampling

As described in sec. 1.4 the central problem with pulsars in binary systems is the varying time delay introduced by a changing distance of the source. The two-body problem has a well-known analytical solution in celestial mechanics, which is the movement on an ellipse around the system’s center of gravity. The precise motion of the pulse-emitting primary star is described by six Keplerian parameters plus post-Keplerian parameters. Most of these are irrelevant for the signal detection because they hardly affect the received pulse train in our case (sec. 2.4.2). For the UCXB system considered in this search we can even neglect the eccentricity (see sec. 2.4.1) which reduces the number of orbital parameters by two. Furthermore the semimajor axis a and inclination with respect to the line of sight i have the same effect on the movement along the line of sight, so that we are effectively left with three unknown orbital parameters:

• The projected semimajor axis expressed in light seconds: X = a/c · sin i

• The binary orbital period P_orb

3http://www.top500.org/list/2012/11/

• The orbital phase expressed as epoch of the ascending node passage: T₀ Sometimes this is expressed as T₉₀= T₀+ P_orb/4.

Together these form a set of unknown trial parameters. A fourth search parameter is the barycentric pulse frequency f_p which plays a distinct role because it is traversed by the FFT and not as an independent parameter.

The orbital time delay from eq. (3) now takes the form trec= t_em+ X · sin

 2π

P_orb · (t_em− T₀)



. (4)

What’s actually measured is the arrival time t_rec of each photon, so eq. (4) is to be inverted numerically to find trial emission times in the source frame t_emand resample the whole observation accordingly.

Of course the movement of the observer in our solar system has a similar distorting effect. I corrected for this beforehand with the well-known celestial and spacecraft orbital parameters using the standard software tool faxbary from HEASoft⁴. For a sufficiently precise correction of that effect one needs good knowledge of the source position, which is not necessarily fulfilled in X-ray astronomy. We will discuss this issue in sec. 2.4.4.

2.4 Residual timing errors

The parameter space considered for the time domain resampling is the simplest possible to approximately account for all time delaying effects. It leaves aside several other parameters whose effect is present in binary pulsars but only cause a negligible distortion to the signal.

2.4.1 Eccentricity

Clearly the projected movement on a Keplerian ellipse is not well described by a pure sinusoidal modulation as we consider in our search. Instead one would have to add eccentricity and argument of periapsis and repeal the degeneracy between inclination and semimajor axis. This would enlarge the search parameter space significantly and make the search unfeasable even with the computing power of years ahead.

Luckily the eccentricity in all our candidate systems is expected to nearly vanish. As depicted by Lorimer (2008) and references therein most observed eccentricities for close binary systems are low and strongly decrease for decreasing orbital periods. This is explained by circularisation over time by tidal friction (Zahn, 1977). Specifically for P_orb on the order of hours the eccentricities should be well below 10⁻⁵. In practice these low eccentricities are even hard to be noticed so that most detections of X-ray binary millisecond pulsars only report upper limits on the order of 10⁻³ (e.g. Chakrabarty and Morgan, 1998). While it remains possible that a small fraction of binaries may have been disturbed by a third-body encounter recently enough to be not fully circularized, the assumption is safe for most of our target systems.

2.4.2 Post-Keplerian Parameters

Post-Keplerian parameters describe the movement of bodies in a slightly disturbed two-body system which can be due to non-gravitational or third-body gravitational forces or to effects of

4http://heasarc.gsfc.nasa.gov/docs/software/lheasoft/

general relativity. Such effects appear especially in tight and fast systems such as the UCXBs that we look at. Stairs (2003) gives a review what role Post-Keplerian parameters play for binary pulsars.

We are especially interested in effects that describe deviations from a harmonical time delay modulation, can therefore not be absorbed in our search parameter sets and possibly decrease our sensitivity. Many post-Keplerian effects are connected to the eccentricity and can therefore be neglected right away.

The first possibly influential effect is the Shapiro time delay which describes an increased travel time that photons take when passing near a massive object such as the companion star in our binary system. The largest expected delay occurs in an edge-on system with i = 90^◦ and is on the order of 50 µs. This is large enough to be measured precisely, however constitutes just a small fraction of the fastest expected pulsation periods. It therefore poses no problem to our sensitivity that would have to be considered in this initial search.

The other effect is a changing orbital period due to gravitational waves, which goes along with a changing X and T₀ as well. Applying the formula given in Stairs (2003) to the most extreme of our candidates yields | ˙P_orb| < 5 · 10⁻¹¹s/s. During one of our observation data sets which last up to six hours the change in period accumulates to ∆P_orb=Porb˙ · T_obs= 1.0 µs, which is well below the admissible tolerances on the order of seconds that we find in sec. 2.5. The effects on X and T₀ are at a similar order of magnitude, so that effects of post-Keplerian parameters are no concern in this search.

2.4.3 Varying pulsation period

Slow and steady drifts in the barycentric pulsation period are precisely measured for a very large number of pulsars (see figure 1). Most pulsars have spin-down rates below 10⁻¹⁰s/s while most millisecond pulsars are even below 10⁻¹⁷s/s. The phase shift at the fastest expected spin period of P_p= 0.8 ms caused by such a drift during a T_obs = 6 hr observation can be estimated as ∆φ = ¹₂P˙_p/P_p²· T_obs² = 3.6 · 10⁻³. This is only a small fraction of one pulsation period and will therefore not degrade the signal.

2.4.4 Positioning error

The observer which is in our case the RXTE spacecraft is moving in our own solar system and consequently contributes a time delay to each observed photon. The earth’s orbit around the sun has a semimajor axis of 1 AU which is 499 lts, while RXTE’s low-earth orbit semimajor axis measures 23.2 ltms. While the earth’s movement dominates for long timescales on the order of years, it is rather uniform on observation timescales of hours, where the spacecraft orbit has a stronger distorting effect. We correct for all these well-known movements as described in sec. 2.3, but the problem is that the time variations depend on the projected movement on the line of sight and therefore on the direction of the source. If the direction is not known precisely, the movement cannot be compensated correctly and a residual error remains which is malicious in our case.

To give an estimate of this, consider the RXTE movement with r = 23.2 ltms and a celestial source position estimate that is off by a small angle ∆φ. For the worst-case projection the arrival time error is ∆t = r/c · ∆φ on each side of the orbit. The pulsation signal with period

P_p will overlap by half a phase and cancel if 2∆t = P_p/2. Taking the most extreme condition of P_p= 0.8 ms sets a positioning accuracy requirement of

∆φ  30⁰ (5)

in order to keep the signal undistorted. This matches half the FWHM spatial resolution of the PCA instrument and is indeed much larger than the position uncertainties on our sources which therefore does not pose a problem.

2.5 Search grid

The time domain resampling strategy requires to pick trial parameters of X, P_orb and T₀ which are sufficiently close to any possible actual parameter set so the worst case signal distortion is limited. This is achieved by placing a grid of parameter sets over the possible space constrained in sec. 2.6. Two questions are important here: What is the necessary grid density and what pattern can be used? Denser grids require more computing power, but sparser grids reduce the average and worst-case search sensitivity.

For a reduced search space with near-exactly known P_orb, Dib et al. (2005) give the solution in the two-dimensional X-T₀ space: First they re-parametrise the space with one complex X^∗ value which encodes both phase and amplitude of the orbit. Then the optimal grid is a circle-packing array with the values of X^∗ placed on the corners of equilateral triangles. To determine the grid distance, one has the find how the maximum Fourier power in the signal drops if the parameters are off by ∆X = |X_actual^∗ − X_trial^∗ |. In the relevant regime close to ∆X = 0 the relative power S follows a central peak shaped as

S(∆X) = J₀(∆X · 2πf_pulse)² (6)

where J₀ is the 0^th order Bessel function (Dib et al., 2005). I also verified this behaviour with simulated data. Now, defining the standard X-grid spacing as δX := 1/(2πf_pulse), I can expand the Bessel function for small values of ∆X:

S(∆X) = J₀

∆X δX

2

= 1 −1 2

∆X δX

2

+ O

∆X δX

4!

(7)

The same behaviour follows for ∆T₀ if we derive the corresponding grid size from the complex representation as δT₀:= P_orb/(2π) · δX/X = Porb/(4π²fpulseX):

S(∆T₀) = J₀

∆T₀ δT₀

2

= 1 −1 2

∆T₀ δT₀

2

+ O

∆T₀ δT₀

4!

(8)

We need to extent the grid to the P_orb dimension as there are no sufficiently precise estimates on P_orb of most of our candidate UCXB systems. In fact Dib et al. (2005) gives a rough estimation of the required P_orb grid size δP_orb and so does Middleditch (2004), but both of them did not get the exact prefactor to match my simulations, so I am investigating it here briefly: Let us consider the case of X and T₀ being correct in the center of the observation. Then the error δP_orb is effectively a varying error in T₀ as ∆T₀(t) = (t − t_center)∆P_orb/P_orb. I find the relative

power S(P_orb) by integrating this over the whole observation, using eq. (8) to the second order:

S(∆P_orb) = 1 T

Z _{T /2}

−T /2

1 −1 2

∆T₀(t) δT0

2

dt = 1 T

Z _{T /2}

−T /2

1 −1 2

∆P_orb/P_orb· t δT0

2

dt = (9)

= 1 −1 2

 ∆P_orbT 2√

3P_orbδT₀

2

= 1 −! 1 2

∆P_orb δPorb

2

, (10)

yielding

δP_orb = 2√

3P_orbδT₀

T = 2√

3P_orb T

P_orb

4π²fpulseX . (11)

A comparison of these theoretical values with the signal drop in simulated data is shown in figure 4. Now I collect the results and establish the standard grid size which corresponds to a 50% signal drop from the peak in each parameter respectively:

δX = 1

2πf_pulse (12)

δT₀ = 1

2πf_pulse P_orb

2π 1

X (13)

δP_orb = 1 2πf_pulse

P_orb 2π

1 X

P_orb T 2√

3 (14)

Since the signal S drops quadratically, the worst-case drop in each of these variable which occurs when the real value is right in the center between two trial values amounts to S_min,1D= (¹₂)² · 50% = 12.5%. The required but unknown f_pulse in eq. (12)–(14) is set to a probable maximum value of f_ref = 700 Hz, up to which the calculated sensitivity drop applies. Faster pulsations can still be found but with continuously decreasing sensitivity.

Now we need to find a 3D grid pattern which fills the parameter space efficiently. The difficulty of that arises with δP_orb which depends on P_orb itself and gives the space an intrinsic curvature (see e.g. discussion in Allen et al., 2013, and references therein) and also causes degeneracies between the parameters. This makes the construction of analytic grids difficult enough that stochastic solutions are often preferred. In our case however we can take advantage of the fact that P_orb is approximately known and only extends over a narrow region. Therefore we can choose the parameters to be orthogonal for P_orb in the center and they will remain near orthogonal in the whole search space. For our typical parameter space dimensions I verified empirically that the metric remains indeed near-euclidean.

My finally chosen geometry is the optimum hexagonal grid in the 2D X^∗ space and on each point therein I place a 1D P_orb grid which exactly fulfils eq. (14) with the X and T₀ derived from the respective X^∗. Now neighbours in the hexagonal X^∗ grid are shifted in P_orb direction by varying amounts. The worst case distance of a point to his nearest grid point that can occur in this grid is in fact when neighbouring X^∗ values are aligned in P_orb direction. Then the distance is composed of one triangle edge-to-center connection and the Euclidean distance in the P_orb direction, which yields a relative distance d =

q (√

3/3)²+ (1/2)² = ^p7/12 = 0.76.

That corresponds to a maximum signal drop of S_max= ¹₂ · 0.76²= 29% and an average of about Savg≈ 12.5%. If the grid is scaled by a factor g, the average sensitivity will scale as

S(g) = 1 − 12.5% · g² (15)

for small values of g. 2.

1185 1190 1195 1200 1205 1210 1215 P

[s]

0.0050 0.0055 0.0060 0.0065 0.0070

X [s]

simulation data 10000s_SP1.4_OP1200_X0.006_pF0.05_I200_BKG0.FITS

0 150 300 450 600 750 900 1050

P (Leahy)

Figure 4: Fourier power (Leahy normalized) at each point on the search grid around the correct parameter set for a simulated typical observation with 200 events/s over 10000 s with 5% pulsed fraction.

This plot shows a 2D cut along the P_orb and X dimension, while T₀ is held at the correct value.

The variation of δP_orb along X is visible. The big ellipse shows the range of theoretical 50%

sensitivity loss according to eq. (12) and eq. (14). The values are actually a bit higher which is due to the second order approximation in eq. (7) and additional random noise. Further from the correct parameter set the power decays less than the central peak approximation used here.

The next step is to determine on optimum grid scaling as a ratio of the identified standard grid size. This is a matter of computing power vs. maximum signal degradation. A similar role is played by the arbitrary sampling time for the FFT in sec. 2.7. Grid size and sampling time must be adjusted together to maximise the sensitivity at a given computing power.

We notice that the binning of time samples at short time steps t_s acts as a low-pass filter with a rectangular window profile. The effect on the power spectrum is obtained with the convolution theorem as a multiplication with a wide sinc² function which damps high frequencies as

S(f ) =

sin(π f t_s) π f ts

2

, (16)

noting that the damping takes effect even below the Nyquist frequency 0.5/t_s (cf. Van der Klis, 1988). The total damping is a multiplication of eq. (15) and eq. (16). The computing effort on the other hand scales as cpu ∝ g⁻³· t⁻¹_s log(T /t_s) where the log factor is insignificant for the scaling behaviour. Both relative damping S and relative computing effort are shown in figure 5 where an optimum ratio can be identified. A value of S close to one with still reasonable computing effort is chosen with a relative gridscaling g = 1 and a sampling timestep of t_s= 1/4096 s, which is also chosen as a multiple of the typical RXTE PCA binning time to mitigate aliasing effects.

0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 relative gridsize

0.0 0.1 0.2 0.3 0.4 0.5

sampling time / pulse period

0.3 1

3 10

30 100

300

1000

0.3 0.2

0.4

0.5 0.6

0.7 0.8

0.9 0.98 0.95

signal degradation cpu effort

ideal combination chosen

Figure 5: Signal degradation S (continuous contours) and necessary computing power (dashed) as a function of sampling time ts and the grid spacing factor g. The necessary computing power rises steeply as one aims to bring S close to 1. An optimum ratio is found on the dotted line and the final choice is marked with a circle.

For the selection of suitable observation data sets we note the total scaling behaviour of the computing effort cpu ∝ t_FFT· N_sets ∝

∼ T f^s· N_sets by deriving the number of gridpoints from eq. (12)–(14):

cpu ∝ ∆P_orb(X_max³ − X_min³ )f_ref³ f_sT²

g³P_orb³ (17)

2.6 Parameter constraints

The computing effort to probe the parameter space of f_pulse, P_orb, X and T₀ depends directly on its size. In this section we put tight constraints on all of these.

The pulse frequency f_pulse dimension has to be probed on the whole interval from zero up to a maximum value by the nature of the FFT. Precise limits on the maximum expected f_pulse are important because it strongly influences the computing effort by determining the data array size and the necessary grid density in sec. 2.5. The typical range for f_pulse is shown by previous detections in table 1 as 182 Hz–550 Hz while the fastest known pulsar spins at 716 Hz (Hessels et al., 2006). These values may be selection biased towards lower frequencies. However they are not far from the ultimate theoretical limit f_K ≈ 1.25 KHz above which mass would be shed from the equator (Andersson et al., 1999). By setting our Nyquist frequency to 2048 Hz and optimizing the search grid in sec. 2.5 to 700 Hz we are most sensitive to the more likely pulse frequencies while all theoretical values are covered in this search.

The orbital period P_orb is completely unknown for about half of all LMXB candidates in the Liu catalogue (Liu et al., 2007) and without that constraint coherent searches remain infeasible

even with the high computing effort we invest. Many LMXB candidates however exhibit X-ray or optical lightcurve modulations due to the viewing angle. 101 such systems are listed in the Ritter catalogue with varying P_orb accuracies, which often cut down the search space sufficiently. We chose to search P_orb within 3σ limits of the given estimate to have a high probability of including the actual P_orb and account for possible small P_orb variations between different observations.

The projected semimajor axis X = a_p/c sin i has no lower constraint except when i is well known. But as the required computing effort increases with X³ in eq. (17), the upper limit is the main computational concern and one can set the lower limit to zero.

To estimate a_p we can take Kepler’s third law (eq. 1) if P_orb and both masses are known. The primary NS mass is very similar in all systems (e.g. M_p = 1.57 ± 0.35M for MSPs in Zhang et al., 2011). The secondary mass is a priori only constrained to . 1M, but can be estimated as follows: The companion must be filling its Roche lobe to produce enough mass transfer for the high X-ray luminosity. The Roche lobe size R_c then equals the companion’s semimajor axis acand is approximated by Paczyński (1973) as

a_c= 0.46a M_p Mp+ M_c

!1/3

. (18)

Combining this with Kepler’s third law (eq. 1) yields

Porb = 8.9 hr Rc

R

!3/2

Mc

M

!−1/2

(19)

(Verbunt, 1993) which determines P_orb as soon as the companion mass-radius relation r_c(M_c) is given.

I analyse several possibilities for R_c(M_c) to ensure that our parameter space can cover each case. Here I take general constraints of the stellar equations of state rather than detailed binary evolution models (e.g. Rappaport et al., 1982) so that we remain sensitive to all possible systems. The implications of each function on a_p via eq. (19) is depicted in figure 6. The three possibilities are hydrogen burning main sequence star, helium burning star or degenerate white dwarf (Verbunt, 1993). I use following equations of state:

• If the companion is a hydrogen burning main sequence star, an approximate solution is R_c= R · 0.76(M^c/M)^0.78(Rappaport and Joss, 1984). Such a companion star can exist down to 0.09M which corresponds to P^orb = 1 hr.

• For the case of a helium burning companion I take the simple approximation from Verbunt (1993) R_c= R · 0.2M^c/M. The resulting a^p constraints are relatively large, however

this type can often be ruled out by optical brightness constraints.

• A white dwarf companion can exist down to very small radii decreasing with increasing mass. I take non-analytic simulation data from Deloye and Bildsten (2003). He, C and O white dwarfs are considered, from which we show the He curve in figure 6 which gives the widest constraint on a_p. At larger radii there is also a strong dependence on the core temperature T_c, which depends on the evolution history of the system and may be up to 10⁷K. A hotter white dwarf can hold more mass inside a given radius and will have a larger a_p.

For each candidate UCXB with a P_orb estimate I can now get an upper limit on X by using figure 6. The last parameter to be constrained is T₀. It is in principle determined together with P_orb from a measured brightness variation but with two limitations: Firstly I can only

10

10

10

Porbit [h]

10

10

10

10

10

10

X [light seconds]

constant f

H burning companion He burning companion He white dwarf (different T) Cluster simulation

LMXB pulsars radio + X-ray pulsar radio pulsars

radio + gamma ray pulsars

Figure 6: Known binary pulsars and theoretical constraints in the Porb-X diagram. The dots show pure radio and radio+gamma ray pulsars from the ATNF catalogue, LMXB X-ray pulsars like the ones we are searching for and the only known LMXB pulsar visible in the radio domain. The lines give predictions of a_p and are drawn for i = 90^◦ which gives an upper upper limit in X. Lines of constant mass function f_X (eq. (2)) are shown for M_p = 1.5M and M_c= 1M, 0.1M, 0.01M (upper to lower line). Predictions for hydrogen (Rappaport and Joss, 1984) and helium burining companion stars (Verbunt, 1993) are drawn in the most relevant ranges and for He white dwarfs with core temperatures of Tc = 10⁵K, 10⁶K, 10⁷K (lower to upper line) (Deloye and Bildsten, 2003). It is noticeable that no pulsars are found in orbits below 40 min, while large numbers of such systems are predicted in simulations of globular clusters (Rasio et al., 2000). This may partly be due to a selection bias, because the fast and widely used search techniques don’t perform well for low Porb and high X.

link brightness fluctuation and orbital time delay phases if I make risky assumptions about the binary geometry. Secondly a well known T₀ constraint at one point is lost after some time t by any inaccuracy in P_orb as ∆T = t · ∆P_orb/P_orb. For example at P_orb = 30 min and

∆P_orb/Porb = 0.5% this timespan is a few days, which is much less than the typical spacing between different observations. I could try not to transfer T₀ into our observation timeframe as assumed in sec. 2.5. But then T₀ would correlate with P_orb making a much finer P_orb grid spacing necessary and undoing the benefit of the T₀ constraint. Therefore I can hardly take advantage of a previously known T₀ and have to search a complete T₀ interval at the length of P_orb.

2.7 Fourier analysis

When a set of trial parameters X, P_orband T₀is sufficiently close to the real values, the resampled time sequence resembles the original coherent pulsar signal which has periodic modulations at the pulse frequency f_p.