Further Evidence for the ∼9 s Pulsation in LS 5039 from NuSTAR and ASCA

The present study aims to reinforce the evidence for the ∼9 s pulsation in the gamma-ray binary LS 5039, derived from a Suzaku observation in 2007 and a NuSTAR observation in 2016. Through a reanalysis of the NuSTAR data incorporating the orbital Doppler correction, the 9.0538 s pulsation was confirmed successfully even in the 3–10 keV range, where it was undetectable previously. This was attained by perceiving an energy-dependent drift in the pulse phase below 10 keV and correcting the pulse timing of individual photons for that effect. Similarly, an archival 0.7–12 keV data set of LS 5039, taken with the ASCA Gas Imaging Spectrometer in 1999 October, was analyzed. The data showed possible periodicity at about 8.882 s, but again the energy-dependent phase drift was noticed below 10 keV. By correcting for this effect, and for the orbital Doppler delays in the LS 5039 system, the 2.8–12 keV periodicity became statistically significant at 8.891 ± 0.001 s. The periods measured with ASCA, Suzaku, and NuSTAR approximately follow an average period derivative of Ṗ≈3.0×10−10 s s−1. These results provide further evidence for the pulsation in this object and strengthen the scenario proposed by Yoneda et al. that the compact object in LS 5039 is a strongly magnetized neutron star.


INTRODUCTION
Through hard X-ray observations with Suzaku in 2007 and NuSTAR in 2016, Yoneda et al. (2020), hereafter Paper I, derived evidence for a ∼ 9 s pulsation from LS 5039, the prototypical gamma-ray binary.This apparently settles the controversy (e.g., Dubus 2013) about the nature of the compact object in this system, indicating that it is a neutron star (NS), rather than a black hole.The result further led Yoneda et al. (2020Yoneda et al. ( , 2021) ) to propose that this NS is in fact has a magnetarmaxima@phys.s.u-tokyo.ac.jp class strong magnetic fields, which are responsible for the particle acceleration and the consequent bright MeV gamma-ray emission from LS 5039.It is also suggested that magnetars can reside, not only as isolated objects, but also in binaries.
In spite of these rich implications, Paper I left several basic issues unsolved.(i) The pulsation in the NuS-TAR data may not have solid significance yet.In fact, Volkov et al. (2021) reconfirmed the Suzaku pulses at the same period as in Paper I (see Table 1), but failed to confirm the pulsation in the present NuSTAR data.(ii) For some unknown reasons, the pulses become invisible in the NuSTAR data below 10 keV, even though no particular feature is seen in the X-ray continuum at ∼ 10 keV (Takahashi et al. 2009;Kishishita et. al. 2009;Volkov et al. 2021).(The Suzaku HXD data are limited to > 10 keV.)This problem is indeed puzzling, and is invoked as an argument against the reality of the ∼ 9 s pulsation (Volkov et al. 2021).(iii) The orbital solutions derived with Suzaku and NuSTAR are not yet fully consistent with each other, particularly in the eccentricity e and the X-ray semi-major axis a x sin i. (iv) The 10-30 keV pulse fraction is much higher in the HXD data, 0.68 ± 0.14, than in the NuSTAR data, 0.135 ± 0.043.
In the present paper, we aim at solving (i) and (ii) above.They are to some extent coupled; any hint of the same pulsation, if confirmed in softer energies, will enhance the pulse reliability, and open a way to the pulse searches utilizing archival soft X-ray data of LS 5039.As to (ii), understanding the origin of the pulse disappearance below 10 keV would provide important information as to the X-ray emission mechanism of this object.
For the above purposes, we first reanalyze, in § 2, the same NuSTAR data as used in Paper I and Volkov et al. (2021), but focusing on energies below 10 keV.Based on that result, we analyze in § 3 the 0.7-12 keV data of LS 5039 acquired in 1999 with the ASCA GIS, hoping to reconfirm the pulsation.In these studies, searches for periodic signals are carried out utilizing the standard epoch-folding method.The epoch-folded profiles are examined for the periodicity significance, employing the Z 2 m statistics (Brazier 1994;Makishima et al. 2016) with the harmonic number set to either m = 2 or m = 4 (Paper I).Since Z 2 m is less widely known than the more standard chi-square technique, in Appendix A we compare the two statistics, and show that Z 2 m much more suited to the present aims.The orbital period of LS 5039 is fixed at 3.90608 d (Aragona et al.2009).All errors refer to 68% confidence limits.

A brief review of Paper I
To begin with, let us briefly review the results in Paper I.Although the observation and basic data reduction described therein are not repeated here, some important numbers are reproduced in Table 1.There, we present not only the Z 2 4 values, but also the reduced chi-square values, χ r , calculated with 20 bins (ν = 19 degrees of freedom) and 16 bins (ν = 15).
The Z 2 4 periodograms (PGs) in the 10-30 keV energy range, from Suzaku and NuSTAR, are reproduced in Figures 1(a) and 1(b), respectively.Like in Figure 3 of Paper I, they are Doppler corrected, using the orbital solutions obtained in Paper I. However, unlike in Paper I, at each period we readjusted a x sin i, the argument of perigee ω, and the initial orbital phase τ 0 (zero when the X-ray object is at the periastron), instead of fixing them to the optimum single values.These parameters were varied over their respect error ranges (Paper I), while e is fixed at 0.306 for (a) and 0.278 for (b).The PG peaks represent the pulse periods on these occasions (Table 1).For the NuSTAR observation, it is given as P NS = 9.05381(3) s. (1) Due to the parameter readjustment, the peaks in  The ordinate is logarithmic, to make the pulse fraction directly comparable.A yellow arrow indicates the pulse peak/bottom ratio in the 10-30 keV profile (Figure 1d).(b) Photons detected with NuSTAR in 3.0-22.1 keV, sorted into two dimensions.Abscissa is the same pulse phase as in (a), while ordinate is the energy in logarithmic intervals which increase by a factor of 1.33.The color coding is black-blue-purple-red-orangeyellow, from lower to higher counts.
files after applying a running average, with weights of 1/4, 1/2, and 1/4 for consecutive three bins.As a result, the statistical error in each bin becomes 0.61 times the Poissonian fluctuation (Makisima et al. 2021a).

A glance at the data below 10 keV
We proceed to the NuSTAR data analysis in soft Xray energies below 10 keV.Starting from the fiducial 10-30 keV PG in Figure 1(b), the lower energy boundary E L was gradually decreased, and the PG calculation was repeated.Then, the maximum Z 2 4 at Equation (1) quickly decreased, from 66.9 in Figure 1(b), to 51.5 at E L = 8 keV, and 28.7 at E L = 6 keV.Likewise, no major peaks exceeding Z 2 4 ≈ 22 were seen in the 7-10 keV or 5-7 keV PGs, even when expanding the period range to P = 9.05 − 9.06 s.We thus reconfirm a conclusion of Paper I, that the pulsation disappears below 10 keV of the NuSTAR data.As described there, a typical pulse fraction is ≈ 0.14 and < 0.03, in energies above and below 10 keV, respectively.
To see what is taking place in softer energies, we next folded the NuSTAR data in several energy bands below 10 keV, employing the same orbital solution as used to produce Figure 1(d).The results given in Figure 2(a) visualize the pulse suppression in lower energies; the pulse amplitudes, if any, are much lower than that in 10-30 keV (vertical yellow arrow).Nevertheless, the charac-teristic "three-peak" structure revealed in Figure 1(d) still remains, with reduced amplitudes, in the 8-10 keV and 6.5-8 keV profiles.Furthermore, from 6.5-8 keV to 8-10 keV, the peaks appear to shift toward later pulse phase.Similar "hard-lag" behavior is also seen between the 3.3-4.3keV and 4.3-5.3keV profiles.
For a further confirmation, Figure 2(b) accumulates 3.0-22.1 keV photons in two dimensions; the pulse phase (abscissa) and the energy (ordinate), where the n-th (n = 1, 2, ., 7) energy bin has E L = 3.0 × (1.33) n−1 keV.To make the result easier to grasp, each row of the plot is rescaled to have a mean of 0 and a standard deviation of 0.5, so the pulse-fraction information is lost.Above ∼ 10 keV, we find a clear pulse ridge running vertically, accompanied by the two sub peaks at about ±0.25 cycles off.These features reconfirm the pulse profile in Figure 1(d).Below 10 keV, in contrast, red inclined stripes are seen to run toward lower left, in agreement with the suggestions of Figure 2(a).
Based on these features of Figure 2, we infer that the pulsation in LS 5039 is present in < 10 keV as well, but its epoch folding is hampered by an energy-dependent systematic shift in the pulse phase.Indeed, from 9 keV to 3 keV, the pulse phase appears to tour almost a complete cycle.Hereafter, we call this phenomenon "Pulse-Phase Drift", or PPD for short, and regard it as a potential cause of the pulse non-detection below 10 keV.

Formalism
To quantify the PPD effect, we modify the folding analysis; below an assumed "break energy" E b , the pulse phase ψ (0 ≤ ψ < 1) is linearly shifted as where E is the photon energy in keV, and R ppd is a coefficient in units of deg keV −1 .The pulse phase at E > E b is unchanged.A positive value of R ppd means corrections for a "hard lag", like in the present case, wherein a lower-energy photon is given a larger phase delay.Figure 2(b) provides an initial guess of E b ∼ 10 keV, and R ppd ∼ 360 • /(6 keV) =+ 60 deg keV −1 .We first attempted to constrain R ppd .To decouple it from E b , which is likely to be above 9.4 keV as in Figure 2(b), we chose the 3-9 keV interval and performed the epoch-folding analysis, incorporating the correction by Equation (2) in which E b is tentatively fixed at 10.0 keV.Scanning R ppd between −90 and +90 (1.5 times wider than the above estimate) with a step of 1.0 (all in units of deg keV −1 ), we studied how the PG peak at P = P NS evolves.In each step of R ppd , we scanned  τ 0 from 0.723 to 0.736 (step 0.001), a x sin i from 47.5 to 48.5 lt-s (step 0.2), and ω from 53 • .5 to 59 4 ≈ +11 means a decrease in the chance occurrence probability by two orders of magnitude.Thus, the correction by Equation (2) appears to be working indeed, yielding two candidates of R ppd .
To estimate E b , and decide between the two R ppd candidates, the same analysis was repeated with the energy range expanded to 3-30 keV.By testing several values of E b from 9 keV to 11 keV, we obtained Figure 3(b).The positive-R ppd peak, which is now at R ppd = 64 ± 2 deg keV −1 , increased markedly, whereas the other candidate diminished.The mechanism working here is instructive.The PPD corrections with R ppd ≈ −15 and ≈ 62 both successfully rectified the < 10 keV pulse phases.Compared to the 10-30 keV pulse template, the soft X-ray profile derived with R ppd = 64 is in a relatively good phase alignment, but that using R ppd = −15 failed to meet this condition.As a result, the positive-R ppd candidate grew up whereas the other became weaker, when we include the 10-30 keV photons which themselves should be insensitive to R ppd .Further consider-

Results of the PPD correction
When the PPD correction as determined above is applied to the photon arrival phases, and the orbital parameters are readjusted only slightly (Table 1), Figure 2 changed into Figure 4.In panel (a), the pulse profiles have become richer in fine structures, and they no longer drift with energy.Likewise, in panel (b), the pulse ridges run mostly straight throughout the broad energy band.The 3-10 keV pulse fraction increased to 0.042 ± 0.022, though still much lower than in 10-30 keV (Figure 2d), and the profiles remain somewhat different from the 10-30 keV one.These are either intrinsic to the source, or due to an inaccuracy of Equation (2) in modeling the PPD effect.
To confirm the pulse recovery in soft X-rays, we created a 3-9 keV PG in Figure 5   peak still remains lower than in Figure 1(b).This is because the 3-10 keV interval, where the pulse fraction is intrinsically lower, contains twice as many photons as the 10-30 keV range, and because PPD-corrected 3-10 profile (black in Figure 4a) somewhat differs from that in 10-30 keV (Figure 1d).
The 3-30 keV PG with no PPD correction, shown in black in Figure 5(c), reveals no peaks at ∼ P NS , because of the dominance of soft X-ray photons.This naturally explains how the NuSTAR pulsation escaped the reconfirmation by Volkov et al. (2021), who started the timing analysis in the 3-20 keV energy interval.

Statistical significance of the PPD effect
When calculating Figure 5(a), we fixed E b and R ppd to the values optimized at P NS .Therefore, this PG is biased toward P NS , and does not correctly reflect the enlarged parameter space.However, if R ppd is also allowed to float at each P , all periods examined will become equivalent, and fluctuations in Z 2 4 at P = P NS , which must be now larger, will provide a reference with which we can evaluate the statistical significance of the peak at P = P NS .(Since E b = 10 keV is outside the energy range, we can fix it.)In Appendix B, we carried out this attempt, and found that the Z 2 4 increase in 3-9 keV at P NS , through the PPD correction (Figure 5b), has a false chance probability of P ch ∼ 5%.
This P ch in 3-9 keV is rather loose, but it changes when the same evaluation is conducted in the 3-30 keV broadband.We obtain a very low chance probability of P ch ∼ 0.004% (Appendix B), for a value of Z 2 4 ≥ 61.51 to appear at P NS (Figure 5c) via the PPD and orbital corrections, in which the parameters except e are all readjusted at each P .Table 2 summarizes these results, together with those derived later from the ASCA data.Thus, the PPD effect, operating below 10 keV of the NuSTAR data, is considered real (see § 4.3).
For reference, the thin red line in Figure 5(c) is a 3-30 keV PG obtained by allowing R ppd to float at each P .(The period search step of 20 µ s is recovered here.)Thus, the effect of the enlarged parameter space is relatively limited in this energy range, typically by ∆Z 2 4 10, which is much smaller than the systematic increase at P NS , ∆Z 2 4 ≈ 42.

Constraints on the orbital parameters
Returning to the 3-9 keV interval, a support (though not quantitative) to the reality of the PPD effect is provided by Figure 6 and Table 1, where we compare the orbital constraints in 10-30 keV (Paper I) and 3-9 keV, derived without and with the PPD correction, respectively.(Here, the orbital parameters are allowed to vary over wider ranges than in Figure 5, for presentation.)The optimum values of ω (panel b) and a x sin i (panel a), derived in 3-9 keV using E b = 10 keV and R ppd = 62, are seen to depend on τ 0 nearly in the same way as the 10-30 keV solution which is free from the PPD disturbance.The difference in P by ≈ 40 µ s is still within relative errors.Thus, the 10-30 keV photons and the PPD-corrected 3-9 keV photons, which are independent, yield nearly the same orbital constraints; they are hence thought to represent the same phenomenon.
A contrasting case is shown in Figure 6 by thin lines denoted as comparison.They represent a typical false solution in 3-9 keV, which we encountered in the calcu- lation of Appendix B. Characterized by P = 10.5285s and R = 76 deg keV −1 , it gives Z 2 4 = 40.6,and its values of a x sin i and ω agree with the 10-30 keV results, if evaluated at a single point of τ 0 ≈ 0.730.Nevertheless, if regarded as a function of τ 0 , this false solution behaves in much different ways from the fiducial one.
Through the reanalysis of the NuSTAR data, we have thus arrived at a scenario that the object is pulsating even in energies below 10 keV, but the pulse phase suffers the PPD effect expressed by Equation (2).This result reinforces the pulse credibility (see § 4.3), and provides an important step forward in solving the issue (ii).It is also suggested that the search for pulsations in LS 5039 can be carried out using rich archival data below ∼ 10 keV, where previous attempts were hampered presumably by the PPD perturbation which was not noticed then.

ANALYSIS OF THE ASCA GIS DATA
Based on the above prospect, we analyze the 0.7-12 keV data of LS 5039 acquired with ASCA in 1999.The aims are to strengthen the evidence for the ∼ 9 s pulsation, and to examine whether the PPD effect noticed in § 2 is present or not.

Observation
LS 5039 was observed once with ASCA (Tanaka et al. 1994), on 1999 October 4, which is 7.94 yrs and 16.92 yrs before the Suzaku and NuSTAR observations, respectively.The gross exposure (the total data span) is T = 63 ks, or 0.19P orb , whereas the net exposure is about 45% of that.We utilize data from the Gas Imaging Spectrometer (GIS: Ohashi et al. 1996;Makishima et al. 1996), placed at the focal planes of the X-ray Telescope (Serlemitsosn et al. 1995).The GIS covers a 0.7-12 keV energy range with moderate angular and energy resolutions.The time resolution is 61 or 488 µs, depending on the telemetry rate.Kishishita et. al. (2009) used these GIS data for spectral and photometric studies, but not for timing analyses.
The data were processed in a standard way.From the identical pair of focal-plane detectors, GIS2 and GIS3, we extracted on-source events over a circular region of radius 5 ′ centered on the source, and co-added the GIS2 and GIS3 events.The arrival time of each event was converted to the barycentric value.Our timing analysis is performed on these photons, without subtracting the background which amounts to 20% and 40% of the total events at 2-3 keV and 10-12 keV, respectively.
The first photon in the data was recorded at MJD 51455.529.It translates to an initial orbital phase of τ 0 = 0.35 or 0.44 (3) based on the Suzaku or NuSTAR ephemeris (Paper I), respectively.The discrepancy between the two values reflects the issue (iii) in § 1.Alternatively, optical observations suggest τ 0 = 0.30 (Kishishita et. al. 2009).Thus, we must treat τ 0 as having a considerable uncertainty.Regardless of this, the final phase of the GIS data is τ 0 + 0.19.Across the observation, the background-inclusive GIS2 +GIS3 count rate in 1-12 keV gradually increased from ≈ 0.2 to ≈ 0.3 cts s −1 .The X-ray light curve of LS 5039 has a good reproducibility (Kishishita et. al. 2009;Takahashi et al. 2009;Yoneda et al. 2023), and this brightening agrees with what is expected for the estimated orbital phase.
We conducted a coarse spectral study by subtracting an appropriate off-source background.The derived 1-12 keV spectrum is featureless, and can be fit adequately by an absorbed power-law model.The photon index is 1.6 ± 0.1, the absorbing column is N H = (6.4± 1.5) × 10 21 cm −2 , and the absorption-removed 2-10 keV flux is (9.7 ± 0.4) × 10 −12 erg cm −2 s −1 .These are typical of LS 5039, and agree very well with those in Kishishita et. al. (2009).

Timing Analyses of the GIS data
The pulse periods from the Suzaku and NuSTAR observations, given in Table 1, define an average pulseperiod change rate as Ṗ = 3.4 × 10 −10 s s −1 (Paper I).A back extrapolation assuming this Ṗ predicts the period at the ASCA observation as P = 8.871 (8.826 − 8.900) s. (4) Here, the upper and lower bounds assume that the actual Ṗ was 1.5 and 2/3 times the above nominal value, respectively.This period tolerance is considered wide enough.In fact, when the pulse-period evolution is ex- pressed as a quadratic function of time t as the upper and lower bounds assumed in Equation (4) imply period-change time scales of Ṗ / P ≈ ±20 yrs which matches the overall observational span of the present study.In addition, the allowance is wide enough to accommodate glitches, in which P would change, in fraction, by ∼ 10 −5 at most (Yu et al. 2013); this is far smaller than the allowance in Equation (4).When searching the GIS data for the pulsation, the parameter space to be examined is huge (see also Volkov et al. 2021), because of the uncertainty in Equation (4), the Suzaku vs. NuSTAR ambiguity in the orbital solutions, the need to consider the PPD effects, and the choice of energy ranges.At the same time, the search grids must be fine enough, because a periodicity in the signal is often accompanied by many side lobes arising via interferences with the quasi-periodic data gaps.Thus, instead of searching at once the whole parameter space using fine grids (spending unrealistic computational times), we proceed in four "Stages" with progressive sophistication and complexity.We can then narrow down the parameter space in stepwise ways, and filter out false side lobes based on their statistical and systematic behavior.This strategy is the same as employed successfully in Paper I, although the actually employed stages are not the same.

Stage 1: Simple periodograms
As the 1st Stage, we tentatively ignore the orbital Doppler effects, and produce Z 2 2 PGs over a period range of 8.2-9.2 s, which includes, and is some 14 times wider than, the range of Equation ( 4).The trial period is varied with a step of ∆P = 100 µs.We use m = 2, because the pulse profiles would still be smeared by the orbital Doppler effects, and lack sharp features.Keeping the upper energy boundary at 12 keV, we tested four values of E L ; 5.5, 6.0, 6.5, and 7.0 keV.
As in Figure 7, a promising result was derived when using the 6.0-12 keV interval; the PG reveals a clear peak, with Z 2 2 = 25.1, at a period of P = 8.8820 s ≡ P 0 (5) which is within the tolerance of Equation ( 4).Since Z 2 2 obeys a chi-square distribution of 4 degrees of freedom, the chance probability for a value of Z 2 2 ≥ 25.1 to appear in a single trial is P (1) ch = 4.8×10 −5 .Given the data span of T = 63 ks, the examined period range contains about T /8.2 − T /9.2 ≈ 835 independent Fourier waves.The post-trial probability hence becomes P ch ≈ 835 P (1) ch ≈ 4.0%.After multiplying this by four, the number of E L tested, we obtain P ch = 16% which is given in Table 2 together with those from NuSTAR.
The above result, 16%, is relatively secure, but not constraining.As already justified, we may consider only those trials which fall in the uncertainty range of Equation (4).Then, P ch reduces to ≈ 1.2% (Table 2), which is sufficiently small.Thus, P 0 is regarded as a start point of our ASCA timing analysis.
When E L > 6.0 keV is used, the PG peak at P 0 diminishes, presumably because the photons would be too few.The peak also decreases when E L is lowered from 6.0 keV, in spite of an obvious increase in the signal photon number.This suggests that the GIS data are also affected, at lower energies, by the PPD perturbation.

Ṗ
Although we have obtained reasonable evidence for a periodicity that is consistent with those from Suzaku and NuSTAR, we still need to perform two timing corrections before concluding on the pulse detection.One is obviously that for the orbital delays, and the other is that for the PPD effect as suggested above.The former, though well formulated, involves multiple parameters (e, a x sin i, ω, and τ 0 ) each with residual uncertainties (Paper I).The latter, on the other hand, has only two parameters (E b and R ppd ), but the formalism is only poorly established.We also need to decide whether to conduct the two corrections simultaneously or in series, and if the latter, which should be the first.
Fortunately, the orbital phase of Equation ( 3) is such that the line-of-sight velocity of the source varies rather linearly with time t (see § 4.1); that is, the pulse period changes with an approximately constant rate of Ṗ7 = 0.92 ± 0.15, where Ṗ7 represents Ṗ in units of 10 −7 s s −1 .The error includes both the uncertainty of τ 0 in Equation (3), and the changes during the observation.We hence remove first the orbital effects approximately, by considering Ṗ in the epoch-folding procedure, in place of the full orbital corrections.The number of parameters then reduces from four to one, which is a big advantage.We also switch from m = 2 to m = 4 (Paper I) because fine structures will emerge in the folded pulse profiles.
As our Stage 2 attempt, we repeated the PG calculation with the GIS data, this time incorporating Ṗ , to obtain Z 2 4 on a (P, Ṗ ) plane.We scanned P over 8.82-8.90s which is comparable to that of Eaquation (4), with an increment of ∆P = 73.2µs, and Ṗ7 from 0.5 to 1.5 with an increment of ∆P 7 = 0.0313.Like in Stage 1, we tested several values of E L , to compromise the signal statistics and the PPD disturbances.Then, E L = 5.3 keV was found to generally maximize Z 2 4 , rather than the 6.0 keV threshold favored in Stage 1.
The result derived from the 5.3-12 keV interval (total 859 photons) is presented in Figure 8(a).Many yellow straight lines run in parallel, with an inclination of ≈ T /2 against the abscissa.Particularly bright are two of them, which are expressed on the plane as P (s) = 8.8820 − 1 2 T Ṗ = 8.8820 − 0.0032 Ṗ7 ≡ P 1 P (s) = 8.8412 − 1 2 T Ṗ = 8.8412 − 0.0032 Ṗ7 ≡ P 2 . (7) They are identified in panels (b) and (c) of Figure 8, which provide two horizontal cuts across panel (a).Of these, P 1 evidently connects to the P 0 peak in Figure 7.
The orbital correction is thus working at least on the P 1 branch, because it becomes highest at Ṗ ∼ 0.7 as expected, reaching Z 2 4 = 32.1 which exceeds the value (Z 2 4 = 22.3) at Ṗ = 0.The other branch, denoted as P 2 , becomes higher than P 1 at Ṗ ∼ 0.9, whereas it gets weaker toward Ṗ → 0, connecting to a weak counterpart in the inset to Figure 7.
Below, we focus on P 1 and P 2 .Either of them is actually a family of parallel lines, with a separation in P by 2-15 ms.This fine structure is created presumably when an enhanced power in the data, from either the pulsation or noise fluctuations, is split by the observing window of ASCA which is a near-Earth satellite.Its data are subject to data gaps, roughly synchronized with the spacecraft's orbital period P sc ≈ 5.5 ks.Hereafter, we collectively call the family near P 1 "Solution 1" (or SOL-1), and that near P 2 "Solution 2" (or SOL-2).

Stage 3: PPD corrections
Now that the orbital effect is successfully approximated by Ṗ , we proceed to Stage 3, namely, incorporation of the PPD correction using Equation (2).Considering that the correction for NuSTAR worked down to the instrumental lower bound of 3.0 keV, the analysis utilizes E L = 2.8 keV, which is justified later.Then, we specify a value of Ṗ in the range of Equation ( 6 The above procedure was carried out for five values of Ṗ7 , from 0.7 to 1.1 with a step of 0.1.The obtained five PGs are superposed in Figure 9(a).A dominant peak reaching Z 2 4 ≈ 52 has emerged at P ≈ 8.89 s, which favors Ṗ7 = 0.9 (red).For reference, we tentatively expanded the period search range further to 7.0-11.0s, as used in Paper I for the initial study of the NuSTAR in 2.8-12 keV, calculated over a broad period range with ∆P = 100 µs, using the PPD correction in which R ppd is optimized at each P (see text).Black, green, red, blue, and gray assume Ṗ7 = 0.7, 0.8, 0.9, 1.0 and 1.1, respectively.(b), (c), and (d) show results over the limited period rage, using ∆P = 20 µ s, for Ṗ7 = 0.9, 0.75, and 1.05, respectively.(e) and (f) give details of the P ′ 2 peak in (d) and the P ′ 1 peak in (b), respectively, where the blue trace presents the optimum R ppd .(g) A superposition of the P ′ 1 peaks for six values of Ṗ (given in color) which cover the range of Equation ( 6).
data, but the peak at 8.89 s still remained the highest, with the next one having Z 2 4 ≈ 49.1.In the same way as in Stage 2, the period search range is hereafter limited to 8.82-8.90s, which accommodates the dominant peak.We re-generated PGs over this range, with a finer step of ∆P = 20 µs.The results are shown in panels (b), (c), and (d) of Figure 9, which employ Ṗ7 = 0.9, 0.75, and 1.05, respectively.The dominant peak in (a) is reproduced in (b) and (c), at a period of P ≈ 8.886 s, which we hereafter denote as P ′ 1 .Although it differs by ∆P ≈ 6 ms from P 1 itself (P 1 = 8.8791 s for Ṗ7 = 0.9), it clearly belongs to the SOL-1 family, because we find (1/P 1 − 1/P ′ 1 ) −1 = 2.07P sc which indicates interference by the observing window.The secondary peak in (a) becomes highest in (d), and seen at P ≈ 8.855 s ≡ P ′ 2 .This P ′ 2 peak belongs to the SOL-2 family, because (1/P 2 − 1/P ′ 2 ) −1 = 0.84P sc holds, Regions around the P ′ 2 peak in Figure 9(d) and around the P ′ 1 peak in Figure 9(b) are expanded in panels (e) and (f), respectively.There, the optimum value of R ppd is shown in blue.Thus, P ′ 1 and P ′ 2 demand R ppd ≈ −23 and ≈ −63 deg keV −1 , respectively.Although the height of P ′ 2 reaches a maximum of Z 2 4 = 44.7 when the involved parameters are trimmed, it remains considerably lower than that of P ′ 1 .From these examinations, we regard SOL-1 as the more likely solution family, among which P ′ 1 is the best pulse-period candidate under the Stage 3 formalism.The P ′ 1 peak in Figure 9(f) is unlikely to be caused by statistical fluctuations, because its chance probably is sufficiently low, P ch ≈ 0.2% (Table 2), as evaluated in Appendix C. We hereafter concentrate on P ′ 1 .The P ′ 1 peak in Figure 9(f) is very sharp, because these PGs are calculated each for a fixed value of Ṗ .To reflect the obvious error propagation from Ṗ to P , Figure 9(g) superposes the P ′ 1 peaks from five values of Ṗ .The composite peak is broader not only than that in Figure 9(f), but also than those in Figure 1 by several times, due to the shorter data span of ASCA.
Figure 9(g) displayed Z 2 4 using P as an independent variable, Ṗ as a secondary parameter, and R ppd as an implicit parameter allowed to vary over a certain range.By interchanging the roles of P and R ppd (but keeping that of Ṗ ), we obtain Figure 10, which is similar to Figure 3. Thus, the 2.8-12 keV GIS data constrain R ppd very tightly.Using these results, and further trimming E b , the maximum pulse significance of Z 2 4 = 53.6 has been obtained under a condition of P ′ 1 = 8.8860 ± 0.0003 s (8a) Ṗ7 = 0.88 ± 0.10 (8b) criteria; it is statistically significant ( § 3.2.1), it lies close to the Suzaku plus ASCA extrapolation, it is consistent with the orbital Doppler effects, and it shares the same E b with the NuSTAR data.We hence regard P fin as the pulse period at the ASCA observation, which is unaffected by the Suzaku vs. NuSTAR ambiguity in the orbital solution.

Additional remarks
To complement the Stage 4 analysis, remarks may be added.First, the peaks in Figure 11(a) and (c) are considerably broader than the composite P ′ 1 peak in Figure 9(g), and by an order of magnitude than the Suzaku and NuSTAR results (Figure 1, Figure 5).This is because P couples strongly with τ 0 (dashed green trace in Figure 11a) under the very limited orbital coverage, and τ 0 is allowed to vary beyond the interval where the Ṗ approximation works.
Next is how the orbital correction works on SOL-2.This is shown in Figure 11  Case 2 : Case 1 plus the orbital correction (Figure 12a).Case 3 : Case 2 plus the PPD correction (Figure 12b).
s, or P ′ 2 + 3 ms, and this increment is consistent with the expected orbital effects.However, the peak, with Z 2 4 = 44.3, is still significantly lower than those in paneld (a) and (c), Z 2 4 ≈ 50.We reconfirm that P fin of Equation ( 9) is the best pulse period at the barycenter of LS 5039.The SOL-2 peak may be a side lobe of P fin .
Finally, we recalculated Figure 11(a), allowing R ppd to float at each P .Then, just like the relation between the blue and red traces in Figure 5(c), the baseline of Z 2 4 was enhanced by ∆Z 2 4 ∼ 10, but the region around P fin remained intact.Allowing E b to float gave no larger effects.These results are not displayed.

Pulse profiles with the GIS
In Figure 12, we compare energy-sorted pulse profiles obtained with the GIS, before (panel a) and after (panel b) the PPD correction.They are all corrected for the orbit using the parameters in Table 1, and folded at P fin of Equation ( 9).In Figure 12(a), the characteristic threepeak structure is reconfirmed in the three narrow-band profiles.At the same time, we notice a clear "phase-lag" property which is reminiscent of Figure 2(a), except that the phase shift has the opposite sign.As a result, the broadband (2.8-10 keV) profile becomes rather dull.As in Figure 12(b), the PPD correction brought the narrowband profiles into an excellent mutual phase alignment, and made them similar to one another, all clearly exhibiting the three-peak structure.The only exception is the 1-2.8 keV result, where the pulses were not confirmed ( § 3.2.3)even when applying the orbit plus PPD corrections.(The purple profile for this softest band is shown without any PPD correction.) Table 3 summarizes how the 2.8-12 keV pulse fraction increases with the timing corrections.The values of Z 2 4 , P ch , and χ r are also given.Thus, the combined application of the orbital and PPD corrections enhances the pulse fraction, and is essential in detecting the pulsation in energies below 10 keV.However, the increase in the pulse significance is not obvious here, because the pretrial probability P (1) ch does not consider the trial number which evidently increases toward Case 3. Interestingly, the final broad-band profile (in black) and the 10-30 keV profile with NuSTAR (Figure 1a) look very much alike.In Figure 12(c) we superpose them together, after appropriately scaling the GIS data (see below).The striking agreement between them strengthens that the GIS and NuSTAR data represent the same phenomenon, i.e., the source pulsation.
In Figure 12(c), the pulse-phase origin of the GIS data was shifted by +0.15.This is because the relative pulse phase cannot be determined uniquely between the two data sets separated by 16.9 yr s.We also rescaled the GIS counts at the j-th bin, C j , into C ′ j , as where C is the average, and the factor 4.27 is the ratio between the NuSTAR and GIS averages.The coefficient 0.54, representing "AC component", means that the 10-30 keV pulse fraction with the NuSTAR is about half that of the 2.8-12 keV GIS data.Indeed, the 2.8-12 keV GIS pulse fraction with the full timing corrections, 0.21 (Table 3), is 1.6 times that of the 10-30 NuSTAR profile (0.135±0.043), although it is much lower than the 10-30 keV Suzaku result of 0.68 ± 0.14.This pertains to the issue (iv) raised in § 1, and suggests that the pulse fraction of LS 5039 can vary significantly, in spite of the stable orbital X-ray light curves.The profiles in Figure 12(b) use the e = 0.306 solution in Table 1.However, even when using the e = 0.278 solution, the profile changes little, rather than becoming similar to the 10-30 keV HXD profile in Figure 1(c).When the SOL-2 parameters (Figure 11b) are employed, the three peaks still persist, but the profile changes considerably, with the left sub peak becoming highest.In the present work, we first reanalyzed the NuSTAR data, focusing on energies below 10 keV.We found that the 9.0538 s pulses lurk in the data even below 10 keV, but their detection is hampered by a systematic pulsephase drift from ≈ 10 keV down to at least 3 keV.Applying an energy-dependent correction to the timings of 10 keV photons for this PPD (Pulse-Phase Drift) effect, the pulses were recovered both in 3-9 keV and in the 3-30 keV broadband, where the pulse detection was difficult before the correction.The PPD effect in the NuSTAR data is considered real, because it has P ch = 0.004% ( § 2.3.3,Table 2) in 3-30 keV.
The orbital solution found with the PPD-corrected signals, either in 3-9 keV or 3-30 keV, agrees with that from the hard-band (10-30 keV) NuSTAR data which are free from the PPD disturbance.Therefore, the X-ray pulses at P 0 are considered to originate from LS 5039, both above and below 10 keV.

Results from the ASCA GIS data
We next analyzed the ASCA GIS data of LS 5039 in four Stages, as depicted in Figure 13 superposed on the predicted orbital Doppler curves.In Stage 1 ( § 3.2.1), a simple m = 2 PG in 6-12 keV, covering a wide period range of P = 8.2 − 9.2 s, gave evidence for periodicity at P 0 of Equation (5) (dashed green line in Figure 13).It has P ch ≈ 16% or ≈ 1.2% (Table 2), when evaluated in the period interval of 8.2-9.2 s, or 8.82-8.90s (consistent with Equation 4), respectively.
The Stage 2 analysis ( § 3.2.2),using the 5.3-12 keV photons, m = 4, and the 8.82-8.90s range, mimicked Alternatively, the 3-30 keV result in Table 2 may be utilized, but it cannot be combined with that in 10-30 keV, because the energy intervals overlap.Then, we return to the 9.025-9.065s period range, which contains 165 independent Fourier waves (for the data span of 3.9 d).Multiplying the 3-30 keV estimate of P ch ∼ 0.004% by 165, we obtain P ch ∼ 0.7%, which is comparable to the first estimate.In either case, the probability for the NuSTAR pulse detection to be false diminishes by an order of magnitude from Paper I, with the pulse confidence increasing to > 99%.
We next consider the ASCA data analysis, where the significance of the ∼ 8.9 s pulses was evaluated at Stage 1 and Stage 3. Of them, we choose the Stage 3 estimate, P ch ∼ 0.2%, because it is based on the largest number of GIS photons available for the pulse search.Again, this P ch primarily means a probability for the PPD effect ( § 3.2.3),but it can also be regarded as the probability to observe, through the PPD correction, a value of Z 2 4 > 53.6 by chance.It takes into account the trials in P (over 8.82-8.90s), Ṗ (substituting for a x sin i, ω, and τ 0 ), R ppd , E b , and E L (Appendix C).
Thus, the chance probability of the ∼ 9 s pulsation in LS 5039 is estimated as P ch ∼ 0.15% with Suzaku (Paper I), ∼ 0.35% or ∼ 0.7% with NuSTAR, and ∼ 0.2% with ASCA.We refrain, however, from combining these values of P ch , as they may not be independent of one another.Furthermore, deriving too small values of P ch would be meaningless when considering various systematic uncertainties; e.g., the unsolved issues of (iii) and (iv), the remaining half of (ii), possible biases in Equation ( 4), and the fact that the initial Suzaku search did not cover the period range below 1 s.
In any event, the ∼ 9 s pulsation of LS 5039 is thought to be significant with > 99% confidence in all the three data sets.This affirmatively settles the objective (i), and strengthens the conclusion in Paper I, that the compact object in LS 5039 is a magnetized NS.The mass estimate of 1.23 − 2.35 Solar masses, derived in Paper I assuming the pulsation to be real, is consistent with this conclusion.

Long-term spin down of LS 5039
Taking the pulsation for granted, the pulse periods of LS 5039, measured for 17 yrs from 1999 to 2016 with ASCA, Suzaku, and NuSTAR, are plotted in Figure 14.We derive Ṗ = 2.6 × 10 −10 s s −1 from ASCA to Suzaku, and Ṗ = 3.4 × 10 −10 s s −1 from Suzaku to NuSTAR.Thus, Ṗ is inferred to change mildly, around an average rate of 3.0 × 10 −10 s s −1 which implies a characteristic age of about 480 yr.The system is indicted to be very young, in agreement with the fact that the optical companion is a massive star.Now that not only P but also Ṗ reported in Paper I were reconfirmed, major discussions and conclusions given there remain intact.Namely, the bolometric luminosity ∼ 10 36 erg s −1 of the source can be powered by neither the rotational energy loss nor mass accretion.The stellar winds cannot provide the required energy input, either.Therefore, the compact object in LS 5039 must be a magnetically powered NS, or a magnetar in a binary.The lack of mass accretion from the stellar winds can be explained by assuming that the Alfvén radius exceeds the gravitational wind-capture radius.
The observed Ṗ could be explained by emission of magnetic dipole radiation, if the dipole magnetic field reaches B ∼ 10 15 G (Paper I).However, this cannot be the dominant spin-down mechanism, because it would predict Ṗ to decrease with time, in disagreement with the present result.As a more likely scenario, the object may spin down via interactions with the stellar winds (Paper I), and fluctuations in this process produce the mild variation in Ṗ .At the same time, these interactions presumably cause the magnetar's magnetic energy to somehow dissipate into particle acceleration and the gamma-ray emission (Yoneda et al. 2021).Admittedly, however, this scenario is still subject to many unknowns, including where the X-ray emission in fact comes from, how the NS's magnetic energy is converted to that of the particle acceleration, and how this process is "catalyzed" by the stellar winds.

The nature of the PPD phenomenon
Let us consider possible astrophysical origins of the PPD phenomenon, after Miyasaka et al. (2013) who discussed various possibilities to explain the behavior of GS 0834−430.Except this Be pulsar, LS 5039, and the two magnetars mentioned in § 4.2, this phenomenon is rather rare among compact X-ray sources.Therefore, it may take place under some limited conditions, e.g., the presence of magnetar-class magnetic fields which LS 5039 is suggested to harbor (Paper I).
The PPD observed from LS 5039 has three charcteristics; (1) E b ≈ 10 keV agrees between the two data sets; (2) between them, the sign of R ppd is opposite; (3) the effect suddenly sets in at E = E b , below which the pulse phase depends linearly on E. Of these, (2) rules out reprocessing of harder photons into softer ones after some delays (soft lag), or hardening of softer photons by, e.g., Compton up-scattering (hard lag).In contrast, as proposed by Miyasaka et al. (2013), an energy-dependent beaming of the X-ray emission may work.Then, (1) and (3) suggest that the threshold energy E b has a certain physical meaning, across which the properties of photon emission and/or transfer change distinctly.
One possibility to explain E b is electron cyclotron resonance (e.g., Araya & Harding 1999) in the strong magnetic fields somewhere around the NS.However, as discussed in Paper I in the magnetar context, the stellar winds flowing toward the NS would be interrupted at the Alfvén radius, R A ∼ 2 × 10 10 B 1/3 15 cm, where B 15 is the surface magnetic field of the assumed magnetar in units of 10 15 G.At R A , the magnetic field would decrease to B A ∼ B 15 (R NS /R A ) 3 ∼ 1.3 × 10 7 G regardless of B 15 , where R NS ∼ 10 6 cm is a typical NS radius.This B A is five orders of magnitude lower than would explain the electron cyclotron resonance at ∼ 10 keV.
Putting aside (3), a mechanism to deflect the Xray propagation direction, in an energy-dependent way, is well known in laboratory; namely, X-ray diffraction techniques.Taking a transmission grating as an example, the 1st order diffraction beams are deflected, from the indecent beam direction, by an angle which is inversely proportional to the X-ray energy.This angle has opposite signs between beams of the order +1 and −1, as required by (2).However, in the present astrophysics setting, it is totally unclear what serves as the diffraction grating.Furthermore, the energy dependence in this scenario differs from what is seen in LS 5039.
In this way, satisfactory explanations of the PPD phenomenon are yet to be sought for.Nevertheless, it may possibly give a clue to strong-field physical phenomena that may be taking place in this interesting object.Future X-ray polarimetric observations are encouraged.

The Suzaku versus NuSTAR discrepancy
Although the issues (i) and (ii) have been solved at least partially, we are still left with (iii) and (iv), or in a word, Suzaku vs. NuSTAR discrepancy.As to (iv), at present we can say only that the pulse fraction of this source may vary considerably for unknown reasons.
Let us consider (iii).It is not due to underestimations of errors (Paper I) associated with the orbital parameters.In fact, when e, a x sin i, and ω of the Suzaku solution are forced onto the 10-30 NuSTAR data, the peak Z 2 4 , which was 66.87 (Table 1), worsens to 29.6, even when allowing ample tolerances for P and τ 0 .Vice versa, the peak Z 2 4 of the Suzaku data degrades from 67.97 to 26.4, when e, a x sin i, and ω are fixed to the NuSTAR solution while P and τ 0 are allowed to vary widely.The discrepancy is statistically significant.
The present ASCA results provide a possible clue to this issue.The maximum value of Z 2 4 = 50.7 or 49.9 (Table 1), achieved in Stage 4, is puzzlingly lower than that (Z 2 4 = 53.6)obtained in Stage 3 where the orbital Doppler changes in P are approximated by a constant Ṗ .Obviously, this is opposite to the expectation, that the full orbital correction in Stage 4 should be more accurate, and would give a higher Z 2 4 .Then, the actual Doppler curve might be slightly deviated from those predicted by an ideal elliptical orbit, and happened to be more straight during the ASCA observation.Such perturbations can arise if, e.g., LS 5039 is a triple system with an unseen third body (e.g., Bosch-Ramon 2021), or the X-ray emission region is somewhat (∼ a few lts) displaced from the NS, or the one-day period variation in the radial velocity of the optical companion (Casares et al. 2011) contributes, or the NS undergoes free precession.
Among the above possibilities, we favor the freeprecession scenario, because nearly all magnetars are deformed by their internal magnetic pressure to an asphericity of ǫ ∼ 10 −4 , and hence its rotation and precession periods become different by ≈ ǫ in fraction (e.g., Makishima 2023).The two periods interfere with each other, to modulate the pulse phase at their beat period, T beat ≈ P/ǫ.If ǫ ∼ 0.3 × 10 −4 , we find T beat ∼ 3.4 d, which is comparable to the orbital period (3.906 d) of LS 5039.When this effect is superposed on the orbital modulation, the Doppler curve will vary to some extent from epoch to epoch, and could explain the Suzaku and NuSTAR discrepancies.

CONCLUSIONS
By revisiting the NuSTAR data of LS 5039 after Paper I, and analyzing the ASCA data taken in 1999, we have arrived at the following conclusions.
• In the NuSTAR data at energies below 10 keV, a PPD (Pulse-Phase Drift) effect is noticed with high significance.Its correction recovered the pulsation in 3-9 keV and 3-30 keV, and strengthened the pulse detection with NuSTAR, although astrophysical origin of this phenomenon is unclear.
• When the PPD and orbital corrections are incorporated, the 2.8-12 s ASCA GIS data acquired in 1999 gave evidence for a 8.891 s pulsation.This result, when combined with those from Suzaku and NuSTAR, further reinforces the reality of the ∼ 9 s pulsation in LS 5039.
• For 17 yrs from 1999 to 2016, the object has been spinning down with an average rate of Ṗ = 3.0 × 10 −10 s s −1 , although Ṗ may be mildly variable.
The characteristic age becomes 480 yr.
• Through the new pulse-period measurement and reconfirmation of Ṗ , the present work validates all discussions in Paper I and Yoneda et al. (2021), that the compact object in LS 5039 is a magnetized NS with a pulse period of ∼ 9 s, and is likely to be magnetically powering the non-thermal radiation.
• The pulse-period change along the binary phase could be subject to some perturbations besides the simple Doppler effect in an elliptical orbit.
APPENDIX A: THE CHI-SQUARE AND Z 2 STATISTICS Both Z 2 m (Brazier 1994) and χ r quantify the deviation of a pulse profile from a flat one.While χ r utilizes profiles folded into N bins, Z 2 m operates on discrete photon data, by adding the power from fundamental to the specified maximum harmonic, m.If m = N/2, the two statistics become the same due to the Parseval's theorem.For white noise data, Z 2 m obeys a chi-square distribution with ν = 2m degrees of freedom.
Figure 15(a) compares two PGs from the 10-30 keV NuSTAR data with the orbital correction.The blue one, using Z 2 4 , is similar to Figure 1(b), but the peak is sharper, because we fix the orbital parameters (Table 1).The other PG in brown uses χ r with N = 20 (ν = 19).The two ordinates are rescaled so that their averages and standard deviations become about the same.Although the two PGs look alike, including the peak at P NS , we can point out three advantages of Z 2 m over χ r .First, the peak value of χ r = 4.24 in Figure 15(a) gives a single-trial probability P (1) ch = 1.5 × 10 −9 , which is much worse than that (P  Z 2 4 = 66.87.This is because the power in the pulses of typical X-ray pulsars (except, e.g, the Crab pulsar with very sharp profiles) is limited to the first several harmonics, beyond which the noise power dominates.Hence, χ r is affected by these noise contributions summed up to the Nyquist harmonic.Actually in Figure 15

Figure 1 .
Figure 1.Essentials of Paper I. (a) and (b) are PGs from the 10-30 keV Suzaku and NuSTAR data, respectively.They are corrected for the orbit of LS 5039, with the orbital parameters readjusted at each period (see text).Red/blue shows Z 2 4 with m = 4 (left ordinate), whereas the dashed green curve gives τ0 (right ordinate), with its search range indicated by a green arrow.(c) and (d) are background-inclusive pulse profiles, folded using the orbital solutions specified by the Z 2 4 peak in (a) and (b), respectively.Between (c) and (d), the relative pulse phase is arbitrary.
Figure 1 are broader than in Figure 3 of Paper I, and their widths faithfully represent the uncertainties in P (Paper I).We regard Figure 1(b) as a reference when analyzing the NuSTAR data below 10 keV.Panels (c) and (d) of Figure 1 reproduce the 10-30 keV pulse profiles from Suzaku and NuSTAR, respectively, folded using their best orbital solutions.They are identical to Figure 2 of Paper I, except binning and definition of the pulse-phase origin.Although the two profiles are somewhat dissimilar, they both comprise three peaks per cycle.Here and hereafter, we show pulse pro-

Figure 2 .
Figure2.(a) Pulse profiles of LS 5039 from the NuSTAR data, obtained in 5 energy bands below 10 keV, using the same orbital solution as used in Figure1(d).The ordinate is logarithmic, to make the pulse fraction directly comparable.A yellow arrow indicates the pulse peak/bottom ratio in the 10-30 keV profile (Figure1d).(b) Photons detected with NuSTAR in 3.0-22.1 keV, sorted into two dimensions.Abscissa is the same pulse phase as in (a), while ordinate is the energy in logarithmic intervals which increase by a factor of 1.33.The color coding is black-blue-purple-red-orangeyellow, from lower to higher counts.
a) : The parameters in parentheses are fixed in the analysis.b) : The initial orbital phase τ0 is specific to each observation, and does not need to coincide among the three observations.c) : The Modified Julian Date of the first event in the data, corresponding to τ0. d) : The parameter values of these rows are taken from Paper I. e) : Errors in these rows do not take into account the uncertainties in e.

Figure 3 .
Figure 3. NuSTAR pulse significance in Z 2 4 , shown against R ppd in Equation (2).See text for details.(a) Results in 3-9 keV, where Z 2 4 is in blue and τ0 is in green.(b) The same, but in 3-30 keV, shown for three values of E b ; 10.4 keV (solid red), 10.7 keV (dotted black), and 10.1 keV (cyan).
ing that Figure 2(b) suggests R ppd > 0, and that the value of τ 0 (in green) associated with the positive-R ppd peak in Figure 3(a) is closer to the Paper I result, we adopt this positive-R ppd solution.By trimming E b and repeating the calculation as in Figure 3(b), we obtained an estimate of E b = 10.4 ± 0.3 keV.

Figure 4 .
Figure 4.The same as Figure 2, but the energy-dependent pulse-phase corrections with Equation (2) have been conducted, employing E b = 10.4 keV and R = 63 deg keV −1 .
(a), incorporating the same PPD correction, and the orbital-parameter readjustment at each P like in Figure1(b).Compared to that fiducial PG, the trial period range was expanded by a factor of 13.Then, the highest significance with Z 2 4 = 36.22 is found at P ≈ 9.054 s, whereas other peaks are all Z 2 4 < 34.As detailed in Figure5(b), the peak is indeed right at P NS .As shown therein by an orange trace, the PPD correction also maximizes the Z 2 4 increment at P NS , ∆Z 2 4 ≈ 25.Within ∼ ±1 ms of P NS , this ∆Z 2 4 takes systematically positive values, probably because the pulse power, once restored by the PPD correction, is partially scattered by the observing window.

Figure 5
(c) shows a 3-30 keV broadband PG, produced in a similar way to (b), but using E b = 10.4 keV and R ppd = 64 deg keV −1 .Again at P = P NS , the PG achieves the peak with Z 2 4 = 61.51,which coincides in height with the red peak in Figure3(b).However, the

Figure 5 .
Figure 5. (a) A PPD-corrected 3-9 keV PG with m = 4, using E b = 10 keV and R ppd = 62 deg keV −1 .Like in Figure 1(b), the orbital parameters except e are varied at each P .A vertical green line indicates PNS of Equation (1).(b) Details of (a) around the peak, shown over the same period range as in Figure 1(b).Dashed black line gives the PPD-uncorrected PG, and orange is the difference between blue and black.(c) The same as (b), but in 3-30 keV.The blue trace employs E b = 10.4 keV and R ppd = 64.The red curve shows a result when R ppd is also allowed to float (with E b still fixed) at each P , from −90 to 90 with a step of 1.0.

Figure 6 .
Figure 6.Comparison of the orbital constraints in 10-30 keV (solid lines) and 3-9 keV (dashed lines; E b = 10 keV and R ppd = 62 fixed).The values of P , ax sin i, and ω, which maximize Z 2 4 , are shown as a function of τ0, together with Z 2 4 .The ranges over which these quantities are scanned are indicated by arrows along the ordinates.The lines denoted as comparison are explained in text.(a) Z 2 4 (red) and ax sin i (blue).(b) ω (green) and P (brown).

Figure 7 .
Figure 7.An m = 2 PG from the 6-12 keV GIS data in Stage 1, generated without orbital corrections.The inset shows an expansion leftward of the peak at P = P0 = 8.8820 s.The green arrow represents Equation (4).
), and change P over the broad interval used in Stage 1 (Figure 7) with ∆P = 100 µs.At each P , we scan R ppd with a step of 1.0 deg keV −1 like in the case of NuSTAR, using a wider allowance from −180 to 180 deg keV −1 , and register the maximum Z 2 4 together with the associated R ppd .Tentatively, E b is fixed at 10.0 keV.

Figure 9 .
Figure 9. Stage 3 results from the GIS.(a) PGs with m = 4in 2.8-12 keV, calculated over a broad period range with ∆P = 100 µs, using the PPD correction in which R ppd is optimized at each P (see text).Black, green, red, blue, and gray assume Ṗ7 = 0.7, 0.8, 0.9, 1.0 and 1.1, respectively.(b), (c), and (d) show results over the limited period rage, using ∆P = 20 µ s, for Ṗ7 = 0.9, 0.75, and 1.05, respectively.(e) and (f) give details of the P ′ 2 peak in (d) and the P ′ 1 peak in (b), respectively, where the blue trace presents the optimum R ppd .(g) A superposition of the P ′ 1 peaks for six values of Ṗ (given in color) which cover the range of Equation (6).

Figure 11 .
Figure 11.The final (Stage 4) PGs from the 2.8-12 keV GIS data, in which ax sin i, ω, and τ0 are readjusted at each P , whereas the PPD parameters are fixed (see text).The red and green traces show the maximum Z 2 4 and the optimum τ0, respectively.(a) A result for SOL-1, assuming e = 0.306, E b = 10.0 keV, and R ppd = −23.4.(b) A result for SOL-2, using e = 0.306, E b = 10.5 keV, and R ppd = −62.(c) The same as (a), but employing e = 0.278 (from Suzaku), together with E b = 10.1 keV and R ppd = −23.6.
(b), which was obtained in the same way, but employing E b = 10.5 keV and R ppd = −62 deg keV −1 .The SOL-2 peak now appears at 8.858
Summary of the data analysis 4.1.1.Results from the NuSTAR data

Figure 13 .
Figure 13.A summary of the period estimates with ASCA, superposed on the radial X-ray Doppler curve from the Suzaku (red) and NuSTAR (blue) solutions.Periods and period differences are all shown in units of ms.The periods found in Stage 1 through Stage 4 are indicated by circled numbers.The orbital phase covered with ASCA is indicated in green, and the optimum Ṗ7 by a dotted black line.

Figure 14 .
Figure 14.The pulse period of LS 5039 measured with ASCA, compared with those from Suzaku and NuSTAR.The red line connects the Suzaku and NuSTAR data points, where the blue one represent an average spin-down trend.

Figure 15 .
Figure 15.Comparison of the χr and Z 2 m statistics.(a) The orbit-corrected NuSTAR PGs in 10-30 keV, calculated by Z 2 4 (blue) and χr (brown; ν = 19).(b) The maximum Z 2 m (black) from the 10-30 keV NuSTAR data, shown as a function of m.The associated P (1) ch is given in orange, where the vertical brown arrow indicates typical values obtained with χr.(c) A scatter plot between χr and Z 2 4 , for the results in panel (a), Table 1, and Table 3.The white box with black bars indicates the average and standard deviation for white noise.The values of P (1) ch are shown in orange.

Table 1 .
The orbital and PPD parameters of LS 5039, derived with Suzaku, NuSTAR, and ASCA.a)

Table 2 .
A summary of the statistical significance of the pulsation.These are intermediate periods, and are somewhat different from the final periods given Table1.d) : Not shown because summed values of Z 2 m are used.e) : The orbital Doppler effects are emulated by a constant Ṗ .
a) : The items with † are utilized in the final significance estimation in § 4.3.b) : The method used in the significance evaluation."MC" means Monte-Carlo simulations, and "Control" means a use of the same data over different period ranges.c) :

Table 3 .
Pulse fraction in the 2.8-12 keV GIS data.