1. Introduction
Neutron stars in low-mass X-ray binaries (LMXBs) are among the brightest sources in the X-ray sky. Their emission is predominantly powered by the accretion of matter from the Roche lobe filled late-type companion star; namely, the conversion of the gravitational potential energy of the infalling material into radiation near or on the surface of the compact object. As a result, such neutron stars either emit continuously (persistent sources) or only when the accretion mechanism is reinstated (transient sources). There are plethora of characteristic observational imprints in their X-ray emission, such as twin quasi-periodic oscillations (QPOs) in the kHz regime (van der Klis et al. Reference van der Klis, Wijnands, Horne and Chen1997), thermonuclear X-ray bursts with oscillations (Watts Reference Watts2012). However, the majority of nearly 200 known neutron stars in LMXBs do not emit coherent X-ray pulsations.
Accreting millisecond X-ray pulsars (AMXPs) has emerged as a subclass of neutron star in LMXBs in the last two decades (Patruno & Watts Reference Patruno and Watts2012). These transient systems exhibit coherent pulsations with periods shorter than 
${\sim}10$ ms. The remarkable capability of AMXPs to turn the accretion energy into pulsations raises a question regarding the reason behind the lack of X-ray pulsations from the majority of the LMXBs. There might be no emerging pulsations from majority of LMXBs possibly because of insufficient magnetic field strengths to channel the accreting matter to the magnetic pole. Alternatively, the beamed signal could already be present but weak; it might be further weakened to non-detection while emerging pulsed radiation undergoes various dynamical or physical processes. Detection of intermittent episodes of coherent pulsations in the persistent emission phase from couple of systems (Aql X–1, Casella et al. Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008; SAX J1748.9–2021, Altamirano et al. Reference Altamirano, Casella, Patruno, Wijnands and van der Klis2008) support the case that all neutron stars in LMXBs are likely to radiate X-ray pulse. These observations provided a unique opportunity to better understand the reason of appearance and then disappearance of the pulsed emission.
There are various neutron star atmosphere and surroundings effects that can reduce the pulse amplitude and, in turn, lead to the lack of pulsations. The most notable scenarios for this situation include light bending effect resulting from the extreme gravity of the compact object (Wood, Ftaclas, & Kearney Reference Wood, Ftaclas and Kearney1988; Özel Reference Özel2009), scattering characteristics of the environment surrounding the neutron star (Brainerd & Lamb Reference Brainerd and Lamb1987; Titarchuk, Cui, & Wood Reference Titarchuk, Cui and Wood2002), and the magnetohydrodynamic instabilities in the accretion flow (Kulkarni & Romanova Reference Kulkarni and Romanova2008). Another possible cause is the binary orbital motion, which could smear the already weak signal to non-detection. Such a dynamical effect would introduce significant implications on the pulsed signal especially in tight binary systems. In such systems, it could, in principle, be possible to recover a pulsed signal in the X-ray data if the binary orbital effects are accounted for.
There have been several extensive studies to correct smearing of pulsations due to dynamical effects from binary neutron stars. One approach, the acceleration search provides a partial correction by splitting the data into short segments and assuming that the orbital acceleration of the neutron star to be approximately constant during the segments of the binary orbit which corresponds to the exposure time of the observation (e.g. Middleditch & Kristian Reference Middleditch and Kristian1984; Anderson et al. Reference Anderson, Gorham, Kulkarni, Prince and Wolszczan1990; Wood et al. Reference Wood1991). Later, Ransom et al. (Reference Ransom, Greenhill, Herrnstein, Manchester and Camilo2001), Ransom, Eikenberry, & Middleditch (Reference Ransom, Eikenberry and Middleditch2002) elaborated this technique to obtain template responses in the frequency–frequency derivative (
$f\hbox{-}\dot{f}$) plane. They correlate the Fourier amplitude and phase responses of the real-time data with template responses that are generated with trial acceleration values. Recently, it was extended by assuming jerk (
$\ddot{f}$) to be constant, instead of the orbital acceleration (Andersen & Ransom Reference Andersen and Ransom2018). By doing that they added one more dimension to their parameter space and conducted their search in 
$f\hbox{-}\dot{f}\hbox{-}\ddot{f}$ volume which they called jerk search (Andersen & Ransom Reference Andersen and Ransom2018). These methods are widely employed to search for periodic signal in timing data collected in the radio band.
Another approach to detect smeared pulsations is called semi-coherent search. It was first proposed as a technique to detect weak gravitational wave signals (Messenger Reference Messenger2011) and then applied to search for weak X-ray signals from neutron stars in LMXBs (Messenger & Patruno Reference Messenger and Patruno2015; Patruno, Wette, & Messenger Reference Patruno, Wette and Messenger2018). The focus of this method is to detect weak but continuous pulsations. They attempted to achieve this goal by applying a two-step procedure on the X-ray data. In the first stage, they divide data to segments and process each segment with a bank of templates coherently. These templates are produced such that they account for the Doppler modulation in the phase of the signal by a Taylor expansion in frequency. In the second stage, they incoherently combine the coherent signal power results obtained for each segment. They applied this technique to 12 LMXBs, and no evidence was found for a previously non-detected weak pulsation.
Fundamental characteristics and RXTE observational details of the systems investigated.
a There are two different burst oscillation frequencies reported for EXO 0748–676.
b Average count rates are calculated in the energy range of 
$\sim$3–27 keV from the all available X-ray data excluding thermonuclear bursts, that is, segment from 20 s before till 200 s after the peak of X-ray bursts.
References– 1. Villarreal & Strohmayer (Reference Villarreal and Strohmayer2004); 2. Galloway et al. (Reference Galloway, Lin, Chakrabarty and Hartman2010); 3. Altamirano et al. (Reference Altamirano2010); 4. Markwardt, Strohmayer, & Swank (Reference Markwardt, Strohmayer and Swank1999); 5. Strohmayer et al. (Reference Strohmayer, Zhang, Smale, Day, Swank, Titarchuk and Lee1996); 6. Smith, Morgan, & Bradt (Reference Smith, Morgan and Bradt1997); 7. Muno et al. (Reference Muno, Fox, Morgan and Bildsten2000); 8. Bhattacharyya et al. (Reference Bhattacharyya, Strohmayer, Markwardt and Swank2006); 9. Chevalier & Ilovaisky (Reference Chevalier and Ilovaisky1991); 10. Zhang et al. (Reference Zhang, Jahoda, Kelley, Strohmayer, Swank and Zhang1998); 11. Welsh, Robinson, & Young (Reference Welsh, Robinson and Young2000); 12. Casella et al. (Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008); 13. Wijnands, Strohmayer, & Franco (Reference Wijnands, Strohmayer and Franco2001); 14. Strohmayer et al. (Reference Strohmayer, Zhang, Swank, White and Lapidus1998); 15. Strohmayer & Markwardt (Reference Strohmayer and Markwardt2002); 16. Kaaret et al. (Reference Kaaret2002); 17. Galloway et al. (Reference Galloway, Muno, Hartman, Psaltis and Chakrabarty2008); 18. Thompson et al. (Reference Thompson, Rothschild, Tomsick and Marshall2005); 19. Hartman et al. (Reference Hartman, Chakrabarty, Galloway, Muno, Savov, Mendez, van Straaten and Di Salvo2003); 20. Wachter et al. (Reference Wachter, Hoard, Bailyn, Corbel and Kaaret2002); 21. Kaaret et al. (Reference Kaaret2007);
Physical mechanism behind pulse smearing is the relativistic Doppler effects caused by the rotational motion of the neutron star around the common centre of mass of the binary system. Effectively, Doppler effects modulate the signal by causing delays on the photon arrival times depending on the relative position of the neutron star on its orbit. This delay was formulated by Blandford & Teukolsky (Reference Blandford and Teukolsky1976) for observing and testing various relativistic effects by making use of the information that pulsars are reliable and precise clocks. In this study, we use this formulation to revert the effect of the binary orbital motion to possibly strengthen the pulsed signal. Smearing can redistribute the signal power to other frequencies around the actual frequency such that the power of the smeared signal lies near or below the detection threshold. We, therefore, applied a systematic investigation by first correcting the arrival times of photons including pulsed signals that are just below the detection level, and then searched for periodic signal in the corrected data. Note that earlier strategies to search for weak pulsations relied on sophisticated templates to model the binary effects on the coherent pulsation signal. In our method, we first account for the smearing effect by transforming the photon arrival times to an inertial frame at the common centre of mass of the binary system and then use a simple sinusoidal model to measure the strength and the frequency of the pulsed signal.
Here we present the results of our extensive search for transient episodes of pulsed X-ray emission in the entire Rossi X-ray Timing Explorer (RXTE) data of 13 binary neutron stars systems, as well as our investigations to constrain the orbital characteristics of these systems.
2. Methodology
We performed the search for weak pulsations in two steps. In the first tier, we conducted pulse search on the barycentred data and determined the candidates. In the second step, we applied arrival time corrections to these candidates to account for the effects of the binary motion and researched for pulsations in the corrected data.
2.1. Observations, the 
$\textit{\textbf{1}}{\textbf{\textit{st}}}$ tier search, and results
For the first tier pulsation search, we selected 13 neutron star LMXBL and used all available event mode archival RXTE data. Note that none of these sources show X-ray pulsations above the noise level during their persistent emission phase (except Aql X–1). However, ten of them show prominent QPOs just before, during, or just after they ignite a thermonuclear X-ray burst. Remaining three also have reported periodicities in their burst observations. Nevertheless, these sources (XTE 1739–285, A 1744–361, GS 1826–238) show either only one burst or the burst oscillations are tentative. We assume that thermonuclear burst oscillations correspond to an X-ray emitting hotspot, and the spin frequency of the underlying source is around the frequency of these oscillations (see Watts Reference Watts2012). Binary orbital periods of six of these sources are known either from periodicities or eclipses observed in the X-ray data. We present the list of the sources investigated and their burst oscillation frequencies in Table 1.
Before applying the first tier search, we generated the light curve of each RXTE pointing in the 2–60 keV energy band with 0.125 s time resolution to search for thermonuclear X-ray bursts. We then created good time intervals for each source by excluding the times of identified X-ray bursts. In particular, we excluded the data starting 20 s before the burst peak till 200 s after. This conservative selection excludes any possible contribution from even the relatively longer duration bursts. Note that the source 4U 1728–34 has a type-II X-ray burster (the Rapid Burster) in its RXTE PCA field of view. For this system, we ignored all observations with type-II bursts present and excluded them from our list of good time intervals. We list the total investigated observing time for each source after the exclusion of the burst intervals in Table 1. Finally, we transferred the photon arrival times to the Solar System barycentre to get rid of the relativistic effects of the moving frame of the detector.
In the first tier search, we fixed channel ranges from one observation to the other rather than fixing energy ranges because energy–channel relation of RXTE has changed during its lifetime. This did not create any problem since we have not compared any observation to the other directly. We applied statistical analysis techniques for comparing any result to the other.
To make our search sensitive for the pulsations that are made up of only hard or low-energy X-ray photons, we carried out the first step of our pulsation search in three energy bands; 
$\sim$3–9 keV (absolute channels of 7–24), 
$\sim$9–27 keV (channel range of 25–70), and 
$\sim$3–27 keV (7–70 channels). Here we note that the energy ranges may vary slightly between the various gain epochs as we have kept the channel range fixed. For each energy band, we constructed a 256 s window at the very beginning of each observation and generated a light curve from the photons that are within this window with a 1/2 048 s binning. This bin size corresponds to a maximum frequency of 1 024 Hz in the Fourier domain and it is above all of the reported LMXB burst oscillation frequencies except XTE 1739–285 for which we used 1/4 096 s binning. We constructed the Leahy normalized power density spectrum (Leahy et al. Reference Leahy, Darbro, Elsner, Weisskopf, Sutherland, Kahn and Grindlay1983) from this light curve and calculated the statistical significance of the maximum power between 
$f_{s}\pm 10$ Hz where 
$f_{s}$ is the reported burst oscillation frequency. Note the fact that intermittent weak pulsations are non-stationary by definition. This makes it harder to decide how to deal with the number of trials while converting the signal powers to its corresponding statistical significance values. We decided to calculate statistical significance values by considering each power spectrum on its own and avoid confusion by providing the power levels besides the significance levels. Therefore, first the single trial probability of obtaining the highest Leahy power is found and then joint probability of having the spectrum is calculated with the number of trials (
$N_{\rm trials}$) equal to the number of frequency bins searched, i.e., the number of frequency bins between 
$f_{s}\pm 10$ Hz. At last, significance level of this joint probability is calculated from a normal distribution in the light of the central limit theorem. Here we have calculated the Gaussian significance of the joint probability by considering a two-tailed test considering the p-value. It is to be noted that the significance values do not change by an appreciable amount if one-sided test (considering one tail of the Gaussian) is applied. We then slided the 256 s interval by 16 s, repeated the same procedure to obtain the significance for that time interval, and continued until the end of the observation. This procedure would facilitate the detection and strengthening of any signal which is present of a short time duration.
After calculating the statistical significance of the strongest pulse for each window, we selected candidates by applying the following continuity, coherence, and minimum significance criterion. We require a pulsation candidate for further detailed analysis to have at least 2.5
${\sigma}$ statistical significance for four consecutive time segments and having the maximum Leahy power between 
$f_{s}\pm 2$ Hz. By setting this criterion, we aimed to uncover coherent periodicities around the burst oscillation frequency that are just below the detection threshold. Our first tier search resulted in 75 episodes of pulsation candidates from 10 sources. We also found that the search in the broader energy range resulted in the same pulsation candidates in the lower and upper energy bands. We, therefore, continued our investigations within the absolute channel range of 7–70 (
$\sim$3–27 keV).
Conventional fast Fourier transform (FFT) can only be applied to discrete data with equally spaced time axis. For that reason, photon arrival times are binned and histograms are created to apply FFT. Even though this approach is computationally very effective, it also has downsides. To be able to bin the data one need to choose the starting time of the binning according to the first photon of interest. Power spectrum may slightly change depending on the choice of the starting time of the binning because shifting the starting and ending times of the histogram bins can cause the photons to be redistributed which will change the values of the histogram. For these reasons, we confirmed the pulsation candidates by using a 
$Z^2$ test (Buccheri et al. Reference Buccheri1983) which is computationally more expensive, but it does not require the data to be binned. The 
$Z^2$ power is formulated as
\begin{equation}Z^2_n = \frac{2}{N} \sum_{m=1}^{n} \left[\left\{\sum_{j=1}^{N} \cos (m\phi_j)\right\}^2 +\left\{\sum_{j=1}^{N} \sin (m\phi_j)\right\}^2 \right],\end{equation}
where n is the number of harmonics, N is the total number of photons, and 
$\phi_j$ is the phase of the 
$j{\rm th}$ photon. 
$Z^2$ powers are distributed with a 
${{\chi}^2}$ distribution with 2n degrees of freedom where n corresponds to the desired number of harmonics. In our case, we have chosen the number of harmonics to be 1 since strong harmonics are not reported for the sources of interest. The 
$Z^2$ powers are calculated with a 1/512 Hz frequency resolution between 
$f_{s}\pm 2$ Hz. Significance values are calculated with the same approach as above. We present the results of our first tier search and the corresponding pulsation candidates in Table 2. We present an example case for a pulsation episode in Figure 1. As in this case, none of our candidate pulsation episodes contain burst emission.
Results of the first tier pulsation search.
a Time (UTC) is the starting time of the candidate pulsation episode.
b 
${t_{\rm dur}}$ is measured from the starting time of the first time segment to the ending time of the last time segment that fits our detection criterion of 2.5
${\sigma}$ statistical significance in four consecutive time segments.
c Leahy power, 
${Z^2}$ power, 
${f_{L}}$ and 
${sig_{L}}$ are obtained from the time segment where the strongest pulsation is observed within the corresponding pulsation episode.
d Count rates during the candidate pulsation episodes in the energy range of 
$\sim$3–27 keV.
e This is the pulsation that is shown in Figure 1.
f This is already reported as an intermittent pulsation by Casella et al. (Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008).
Example of a candidate pulsation episode from 4U 1636–536 that lasts for five consecutive time windows. The pulsation candidate and its properties are indicated in boldface in Table 2. (top panel) The light curve of the part of the observation that contains 336 s long candidate pulsation episode. The vertical dashed lines that have the same line style correspond to the starting and ending time of 256 s windows from which power spectra are calculated. (bottom plots) 
${Z^2}$ power density spectra of the five sequential 256 s intervals indicated with vertical lines above. The signal at 579.66 Hz is clearly evident in all plots.
2.2. Binary orbital motion corrected search
Until this point, we searched for coherent pulsations in the data that only the barycentric correction was applied. However, binary orbital motion of the neutron star could smear out an already weak signal and make it undetectable. In our study, we used the appropriate relativistic orbit model (Blandford & Teukolsky Reference Blandford and Teukolsky1976) for the Doppler correction to be able to recover the smeared signal based on plausible orbital parameters. The time delay, 
${t_{d}}$ due to orbital motion is formulated as
\begin{equation} \begin{split} t_{d} = & \left \{ \alpha \left ( \cos E - e \right ) + \left (\beta + \gamma \right)\sin E \right \} \\[6pt]
& \times \left \{ 1 - \frac{2 \pi}{P_{b}} \frac{ \beta \cos E - \alpha \sin E }{ 1- e \cos E } \right \},\end{split}\end{equation}
where 
$\alpha = x \sin w$ and 
$\beta = \left ( 1 - e^{2} \right )^{1/2} x \cos w$ are expressed in terms of the projected semi-major axis, 
$x = a \sin i$, 
${P_{b}}$ is the binary period, E is the eccentric anomaly, e is the eccentricity, w is the longitude of periastron, and 
${\gamma}$ is the term for the gravitational redshift and time dilation. The eccentric anomaly E is defined as
\begin{equation}E - e\sin E = \frac{2\pi t}{P_{b}}\end{equation}
from the Kepler’s equation (Taylor & Weisberg Reference Taylor and Weisberg1989). In our search, we assumed the orbit of the neutron star to be circular which is a plausible approximation for LMXBs (
$e = 0$). This assumption reduced the number of free parameters to three that are the binary orbital period (
${P_{b}}$), projected semi-major axis (x), and the epoch of mean longitude equal to zero (
$T_{0}$). Because of the lack of information and uncertainties about the epoch of mean longitude equal to zero, we decided to apply a 
${\pi/4}$ radian sampling to the phase (
${\psi}$) of the circular orbit that corresponds to eight trial phases in total. This way, we aimed to test roughly all possible configurations that can effect the signal. For the six sources whose binary orbital periods were already reported (see Table 1), we limited the search interval of binary period to within 1 h around their reported value, with a 0.1 h sampling in those 2 h intervals. For the rest of the sources whose binary period is unknown, we chose the range of the trial range of orbital periods to be 2–30.5 h and applied 1.5 h sampling to this time interval. This orbital period interval is expected to correspond to a great majority of LMXB orbital periods (van Haaften et al. Reference van Haaften, Nelemans, Voss, van der Sluys and Toonen2015). Similar to the uncertainties about the epoch of mean longitude equal to zero (
$T_{0}$) of the sources, projected semi-major axis (x) values of LMXB sources were also either unknown or highly uncertain. For this reason, we select a broad range of trial values from 0.01 to 1.91 light-s with 0.1 light-s sampling. Therefore, for each 256 s pulsation candidate time segment, the correction is applied for 
$8\times20\times20$ sets of parameters which span 
${T_{0}}$, 
${P_{b}}$ and 
${a\sin i}$ parameters, respectively. Sampling sizes of the correction parameters were limited mainly due to the computational costs. We tried to cover a wide range of parameters with the smallest grid size possible given the present computational capabilities. We would like to also note that the reported binary periods of our six sources are much greater than 256 s. In practice, one could reduce the parameter space further based on this fact (see (see Ransom et al. Reference Ransom, Greenhill, Herrnstein, Manchester and Camilo2001, Reference Ransom, Eikenberry and Middleditch2002). However, we decided to pursue with the Blandford & Teukolsky (Reference Blandford and Teukolsky1976) correction and applied it in all three dimensions for the sake of completeness and for creating a generalized approach applicable to different sources.
We then searched for pulsations in the binary orbital motion corrected data by employing the Z
$^2$ method between f
$_{s} \pm$ 2 Hz as done before. By doing that, every correction that is applied with a different parameter set gave different statistical significance values for each time segment. Afterwards, we calculated the highest frequency shift that can occur (
${\Delta f_{\max}}$) from the first-order calculations
\begin{equation}\Delta f_{\max} = f_{p} \left | \frac{P_{b}}{{P_{b}}- 2\pi a \sin i} - 1 \right|\end{equation}
and then located recovered pulsations by applying the following set of criterion: We require the recovered pulsations to have at least 3.5, 4.5, and 3.5
${\sigma}$ statistical significance for three consecutive time segments, and we require significance values for these three segments to be higher than before (that is, before the binary motion correction). We choose these significance levels since we wanted to put a significance criterion at least 1
${\sigma}$ higher than the previous (2.5
${\sigma}$) for the detection of the recovered pulses. We require the improvement in the significance levels for three time segments to be achieved concurrently after corrected with the same parameter set. We also require the pulsations to be at the same frequency for these three time segments and the frequency of the recovered pulsation to be at most 
${\Delta f_{\max}}$ away from the frequency of the pulsation before correction.
Results of the binary motion corrected search.
a Time (UTC) is the starting time of the 256 s time window.
b Lines under Aql X–1 are orbital motion corrected results of the previously discovered 150 s pulsation episode by Casella et al. (Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008)
c Detection criterion of recovered pulsations for 4U 1728–34 is different than other sources which is 3.5, 4.0, and 3.5 rather than 3.5, 4.5, and 3.5 for three consecutive time segments. See Section 2.2 for a more detailed discussion.
After applying this set of criterion, we identified recovered pulsations for five sources: EXO 0748–676, 4U 1608–52, 4U 1636–536, Aql X–1, 4U 1728–34 out of 10 that showed pulsation candidates. We obtained degenerate pulsations for almost all time segments that showed recovered pulsations. These pulsations were degenerate such that, for a single time segment there were many recovered pulsations that were obtained with different 
${P_{b}}$, 
${a\sin i}$, and 
${\psi}$ parameter sets. Moreover, pulsation strengths of recovered pulsations for a single time segment were at a similar level. The presence of degenerate recovered pulsations was expected because of the degenerate nature of Doppler effects. For example, to obtain similar Doppler effects for a trial orbital phase there should be a relation between 
${P_{b}}$ and 
${a\sin i}$ such that if 
${P_{b}}$ is short, 
${a\sin i}$ should also be shorter to compensate each other.
Binary orbital period (
${P_{b}}$) of four of the five sources was already known (except 4U 1728–34). This helped us to eliminate the degenerate cases of recovered pulsations. We discarded the recovered pulsations that are obtained with a 
${P_{b}}$ that is different than the reported 
${P_{b}}$. Furthermore, for each time segment we selected the parameter set that gives the highest signal power among the ones that have the appropriate 
${P_{b}}$. Since we did not know the binary period of 4U 1728–34, we adopted a slightly different approach. Since sampling for the 
${P_{b}}$ parameter of 4U 1728–34 is 15 times larger than the ones that the binary periods are known, we decided to decrease the significance criterion of the detection of the recovered pulsation to 3.5, 4.0, and 3.5 from 3.5, 4.5, and 3.5. By doing so, we aimed to report more recovered pulsations and corresponding 
${P_{b}}$, 
${a\sin i}$ couples for such a highly unknown situation. We list the results of the recovered pulsations and the corresponding correction parameters in Table 3. Note that we also report in Table 3 the binary motion corrected results for the previously detected 150 s intermittent pulsation episode from Aql X–1 (Casella et al. Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008). Each line in this table corresponds to a 256 s time window that meets our detection criteria, and it can be seen that some of these time windows are adjacent to each other. We count the adjacent windows as one pulse episode which makes a total of eight new pulsation episodes reported in the table excluding the detection of Casella et al. (Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008).
There were two different spin frequencies in literature that thought to be related to the spin frequency of EXO 0748–676 that are 45 and 552 Hz. We were able to obtain pulsation candidates at both of these frequencies in our first tier search. However, after applying and correcting the photon arrival times, we obtained recovered pulsations only around 45 Hz even though smearing effect of the Doppler modulation is increasing at higher frequencies.
3. Discussion
We have performed one of the most comprehensive search for intermittent pulsation episodes in non-pulsing LMXBs. In particular, we have accounted for the binary orbital motion systematically in our extensive investigations. Recently, Messenger & Patruno (Reference Messenger and Patruno2015) and Patruno et al. (Reference Patruno, Wette and Messenger2018) also accounted for the Doppler shift in searching pulsations from LMXBs with their semi-coherent search. Their work differ from ours at a fundamental level: They aimed to detect continuous but weak pulsations, contrary to our search for intermittent pulsation episodes. Furthermore, eight out of 12 of their LMXB sample did not show burst oscillations. This made their parameter space significantly wider for these LMXBs contrary to more confined spin frequency parameter in our work thanks to the previously reported burst oscillations for the LMXBs. Note that we have nine different LMXBs in addition to four common with their work, namely, Aql X–1, 4U 1636–536, 4U 1608–52, and XTE 1739–285.
Detection of a single 
$\sim$150 s long pulsation episode (Casella et al. Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008) in the entire archival RXTE data of Aql X–1’s persistent emission phase clearly show that the intermittent pulsations are very rare (0.009% of the total 1 645 ks observed time). Note that the pulsation episodes in other sources are also rare with respect to the total observation time: 0.023% for EXO 0748–676, 0.012% for 4U 1608–52, 0.016% for 4U 1636–536, and 0.062% for 4U 1728–34. We should note here that occurrence rate for 4U 1728–34 is expected to be higher since we lowered significance criterion for this source. (see Section 2.2). Besides these recovered pulsations, we also found pulsation candidate episodes that contain weak pulsations around the burst oscillation frequencies for all 13 LMXBs in our sample. These candidate episodes are also short and rare compared to the total searched time of each source: 0.24% for EXO 0748–676, 0.08% for 4U 1608–52, 0.13% for 4U 1636–536, 1.84% for MXB 1658–298, 0.07% for 4U 1702–429, 0.21% for 4U 1728–34, 0.52% for XTE 1739–285, 0.15% for SAX J1750.8–2900, 0.12% for GS 1826–238, and 0.19% for Aql X–1.
Strohmayer et al. (Reference Strohmayer2018) employed the same binary orbital motion correction to account for Doppler smearing in an accreting millisecond X-ray pulsar system, IGR J1762–6143 with an ultra-compact orbit (binary orbital period is 
${\sim}38$ min). They found that the orbital delays weaken the pulsation such that after the correction to its NICER data, 
$Z^2$ power of the pulsation increase up to 
${\sim}196$ which is almost four times the previous value (
${\sim}55$). However, their work differs from ours on a fundamental aspect. IGR J1762–6143 shows strong pulsations as side bands even when no correction was applied. Correcting the photons put the side band structure together and result in a stronger centred pulse. In our approach, we do not observe clear pulsations as side band structures in the non-corrected data since the strength of the pulsations that we have been looking for were not strong enough to be detected after smearing. Nevertheless, it is important to remark that side band structured pulses of IGR J1762–6143 clearly show that smearing effect of the binary orbital motion plays a key role in disappearance of already weak pulsations. This makes it almost mandatory to apply the Doppler correction to the photon arrival times before searching weak pulsations especially for compact (short 
${P_{b}}$) sources considering the high uncertainties in the inclination angles.
We were able to determine orbital inclination of four systems, based upon our criterion that the search in the Doppler corrected data would yield the strongest signal. We suggest the projected semi-major axis (
${a\sin i}$) value of about (5.1–5.4)
$\times10^{8}$ m for EXO 0748–676, and 
$5.7\times10^{8}$ m for 4U 1608–52. For typical orbital separations (a) that can be calculated from the known binary periods, the orbital inclination of these two systems turned out as 29.0–30.9 and 13.9 degrees, respectively. Our search for projected semi-major axis of 4U 1636–536 and Aql X–1 concentrate around two values: 
$1.8\times10^{8}$ or 
$5.1\times10^{8}$ m for the former and 
$3.3\times10^{8}$ or (5.1–5.7)
$\times10^{8}$ for the latter. These correspond to 9.9 or 29.1 degrees of inclination for 4U 1636–536 and 6.2 or 9.6–10.7 degrees for Aql X–1. We would like to note that, donor mass and neutron star mass are assumed to be 0.5 and 1.4 
${\rm M}_{\odot}$, respectively, for calculating the orbital separations. Highly degenerate recovered pulsations and their scattered orbital parameters in the 
${P_{b}}$ – 
${a\sin i}$ plane made it impossible to determine orbital parameter estimates for 4U 1728–34. Two recovered episodes of pulses for EXO 0748–676 were found at the frequencies of 44.5 and 46.4 Hz. We found pulsation candidate episodes around 552 Hz, which is the other suggested frequency, and also there are other studies supporting that (Balman Reference Balman2009; Jain & Paul Reference Jain and Paul2011). However, we did not detect any recovered pulsations after the Doppler correction around the higher frequency. In the light of these results, we argue that the 45 Hz may be the actual spin frequency of EXO 0748–676 rather than 552 Hz.
We have also investigated whether X-ray intensity (that is a measure of mass accretion rate) is different during those pulsation episodes. We compare the average count rate for each source (listed in Table 1) to the source rates during the episodes of coherent pulsations (listed in Table 2). We find no systematic trend in X-ray intensities, neither higher nor lower than their long-term averaged values. We, therefore, suggest that the appearance of pulsations is not linked to any sudden change in the accretion rate.
The lack of pulsations in these bright LMXB sources may not be solely related to the Doppler effects. There are already numerous suggestions attempting to explain the reason of intermittent pulsations. Sporadic behaviour and uncommon nature of such periodicities force the explanations to be somehow related to a rare incident or asymmetry that might be taking place close to the neutron star surface. These explanations can be roughly divided into two categories: one group of explanations assume that the pulsation is temporary while other group assume that the pulsation is always present however so weak that we cannot detect it.
One possible explanation for the lack of pulsations that belongs to the first group is the scattering of the beamed emission by the optically thick media. However, this approach was challenged by Göğüs, Alpar, & Gilfanov (Reference Göğüs, Alpar and Gilfanov2007) where they argued through spectral investigations that the optical depth 
${\tau}$ of the surrounding corona is not thick enough to smear out the pulsations. In order for this explanation to work, temperatures of scattering electrons should be very low (
$\lesssim$10 keV) so that the optical depth would be large enough to screen. The pulsations would then become visible only during a spectral variability, in the presence of a local hole in the screening medium or in the presence of other unique asymmetric geometries (Casella et al. Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008). In such a scenario, one would also expect a dependence of the presence of the pulsation on the accretion rate as in such systems the coronal properties are strong function of the accretion rate. Such a correlation was found to be absent in our analysis.
The second group of explanations is generally related to the magnetic channeling. Neutron stars in LMXBs are weakly magnetized. Therefore, they are typically unable to stop and focus the incoming matter. The appearance of pulsations could either be due to a sudden change in the amount of channeled matter or the strength of their magnetic field could suddenly increase (at least locally) and become capable of channeling matter. The second scenario might occur if there is a sudden decrease in the Ohmic diffusion time which could happen due to a starquake, local disruption of screening currents, or magnetic reconnections (Casella et al. Reference Casella, Altamirano, Patruno, Wijnands and van der Klis2008). It is essential to study multipolar magnetic fields in such neutron stars to find an evidence for strengthening the magnetic fields. Recent serious efforts to understand magnetic field structure of neutron stars by ray-tracing technique (Bilous et al. Reference Bilous2019) may provide new insights on this front.
Another approach within the second group explanations is nuclear burning. We expect nuclear-powered oscillations to last typically a few seconds. However, Strohmayer & Markwardt (Reference Strohmayer and Markwardt2002) reported a long-lasting pulsation during a superburst of 4U 1636–536. This indicates that coherent and longer lasting oscillations may also be related to nuclear burning. This explanation cannot be related to results since we carefully identified the bursts and eliminated their times from our searched sample. There could, in principle, be nuclear burning events which might not be seen radiatively (as thermonuclear bursts), but the stored energy could instead give rise to intermittent coherent oscillations.
Acknowledgements
We are grateful to the anonymous reviewer for conservative comments. We acknowledge support from the Scientific and Technological Research Council of Turkey (TÜBITAK, grant no: 115R034).
Open access
