Modeling the Precession of the Warped Inner Accretion Disk in the Pulsars LMC X-4 and SMC X-1 with NuSTAR and XMM-Newton

The observations used in this analysis consist of two distinct data sets. NuSTAR and XMM-Newton observed LMC X-4 jointly at four different epochs between 2015 October 30 and 2015 November 27. Table 1 lists the observation ID numbers, dates, and exposure times for the LMC X-4 observations. NuSTAR and XMM-Newton also observed SMC X-1 jointly at four epochs between 2016 September 8 and 2016 October 24, and Table 2 contains the information for these observations. Figure 1 shows the one-day averaged MAXI light curves for LMC X-4 and SMC X-1 during the time of observations and indicates the time of observations. For both sources, our observations sample a single superorbital phase.

Zoom In Zoom Out Reset image size

2.1.1. LMC X-4 Data Analysis

We reduced the NuSTAR data for LMC X-4 and SMC X-1 using version 1.8.0 of the NuSTARDAS pipeline and CALDB v20170727. For the XMM-Newton data, we used version 14.0.0 of XMMSAS, with an updated leap-second data file.

For each NuSTAR observation, we used DS9 to select circular source regions with a radius of 120'' centered on the source coordinates. A background region of the same size was selected away from the source. The XMM-Newton observations were taken with EPIC-pn in small window mode to minimize pileup, and we exclusively used the EPIC-pn instrument and not EPIC-MOS for the best timing resolution. In these observations, the source was positioned close to the edge of the chip and we detected small amounts of pileup using the XMMSAS tool epatplot. To mitigate both of these effects, we selected annular source regions with an inner radius set to minimize pileup and the outer radius remaining on the chip. We found typical inner and outer radii for the annular source regions of 13'' and 42'', respectively. We selected XMM-Newton background regions from circular regions of radius 60'', located away from the source area. We filtered all EPIC-pn data to contain only single and double events. We applied a barycentric correction to both the NuSTAR and XMM-Newton data sets using the NuSTARDAS tool barycorr and the XMMSAS tool barycen, respectively. We also corrected the pulse arrival times using the LMC X-4 ephemeris described in Levine et al. (2000).

We show the light curves for the LMC X-4 observations in Figure 2, where the 0.2–12 keV XMM-Newton light curves have been arbitrarily offset from the 3–79 keV NuSTAR light curves for clarity. Several bright accretion flares appear in observations L1 and L4. Brumback et al. (2018) found that these flares contain changes in pulse strength, shape, and phase, which could complicate the relative phase mapping presented in this analysis. For this reason, all flares have been removed from the light curve and only times of direct simultaneous observation between XMM-Newton and NuSTAR have been used.

Zoom In Zoom Out Reset image size

2.1.2. SMC X-1 Data Analysis

We used the same versions of NuSTARDAS and XMMSAS mentioned in the LMC X-4 data analysis to reduce the SMC X-1 data. For each NuSTAR observation, we used DS9 to select circular source regions of 120'' centered on the source coordinates. We selected background regions of the same size away from the source. For the XMM-Newton observations, we used data from the EPIC-pn instrument in timing mode for the best timing resolution. We selected the source region from a column 20 pixels wide, centered on the source. The XMMSAS tool epatplot revealed slight amounts of pileup, and so we excised the brightest central pixel of the source to minimize this effect. We filtered the EPIC-pn data to contain only single and double events. We applied a barycentric correction to both data sets using the NuSTARDAS tool barycorr and the XMMSAS tool barycen, respectively. We also corrected for the pulse arrival times using the SMC X-1 ephemeris described in Falanga et al. (2015).

Figure 3 shows the NuSTAR 3–79 keV and XMM-Newton 0.2–12 keV light curves for the SMC X-1 observations. The light curve from Observation S3 has a very low count rate, indicating that the source was weakly detected. As seen in Figure 1, the superorbital period sampled during our four SMC X-1 observations had a different amplitude and period from the preceding superorbital periods. The variable behavior of SMC X-1's superorbital phase has been previously monitored (e.g., Hu et al. 2011, 2013; Dage et al. 2019) and could be caused by an instability in the accretion disk's radiation-driven warp (Ogilvie & Dubus 2001). The variation in superorbital phase occurring during our observations caused Observation S3 to occur during the low state of the superorbital cycle. We did not detect pulsations during this observation, and we therefore exclude it from further analysis. We also excluded part of the NuSTAR observation for Observation S2, which exhibited changes in pulse behavior not covered simultaneously by XMM-Newton. This change in pulsation is energy-dependent and only affects the NuSTAR pulse profile, so therefore we use the full XMM-Newton observation. The data used in this analysis are plotted in black in Figure 3, and the full NuSTAR observation is shown with reduced opacity.

Zoom In Zoom Out Reset image size

We used epoch folding, via the function epfold found in the Remeis observatory ISISscripts, to find the best period of each LMC X-4 observation, and a Monte Carlo simulation of 500 light curves to find the uncertainties. We used the NuSTAR data in the epoch folding analysis because they were more strongly pulsed than the XMM-Newton data. The best pulse periods are listed in Table 3.

Before creating pulse profiles, we filtered the event files by energy so that the NuSTAR data probed the hard pulsar beam emission (8–60 keV) and the XMM-Newton data captured the soft reprocessed emission (0.5–1 keV). We created energy-resolved pulse profiles using the folding technique in the FTOOL efold, which folds the light curve of each observation by the best period for that observation (see Figure 4). We used 20 bins per phase for these pulse profiles.

Zoom In Zoom Out Reset image size

We found that the Monte Carlo error analysis that we employed for LMC X-4 was not practical for determining uncertainties in the SMC X-1 spin period because of the timing resolution needed to evaluate the ∼0.7 s period. To make the analysis less computationally expensive, we employed the epoch folding technique found in the HENDRICS software (Bachetti 2015) tool folding_search. This epoch folding tool searches the spin frequency and frequency first derivative simultaneously and returns a distribution of ${Z}_{4}^{2}$ statistics (Buccheri et al. 1983). To estimate the 1σ level uncertainty, we fitted this distribution with a two-dimensional Gaussian using the Astropy model Gaussian2D and a Levenberg–Marquardt least-squares fitting routine.

We confirmed that this epoch folding analysis is consistent with the Monte Carlo analysis from LMC X-4 by using folding_search on LMC X-4 observations with high signal-to-noise ratio. We found the results from each method to be consistent, and therefore do not believe that the difference in method will affect our measured pulse periods. The best pulse periods for the SMC X-1 data are listed in Table 4.

Table 4.  Best-fit Spin Periods for SMC X-1 Observations

Observation	Spin Period (ms)	$\dot{P}$ (s s−1)
S1	699.65 ± 0.03	(−1 ± 3) × 10−9
S2	699.59 ± 0.04	(1 ± 3) × 10−10
S4	699.60 ± 0.03	(3 ± 3) × 10−10

Download table as:  ASCIITypeset image

We filtered the SMC X-1 data by energy in the same way as the LMC X-4 data so that our NuSTAR pulse profile captures the hard X-ray component and our XMM-Newton data cover the soft component. We then made pulse profiles with 20 bins per phase (Figure 5) using the Stingray (Huppenkothen et al. 2019) software tool fold_events and the measured period and period derivative.

Zoom In Zoom Out Reset image size

2.3.1. Phase-averaged Spectroscopy

For both data sets, we extracted spectra from the source and background regions described above using appropriate NuSTARDAS and XMMSAS selection tools. However, we did not select background spectra for the XMM-Newton observations of SMC X-1 because the source flux dominates the EPIC-pn timing mode CCD (e.g., Ng et al. 2010).

We grouped all NuSTAR spectra into bins with a signal-to-noise ratio of 18 and all XMM-Newton spectra with a minimum of 100 counts per bin, which produced good statistics. We fitted the phase-average spectra in the range 0.6–50 keV.

We modeled the spectra in Xspec version 12.9.1 (Arnaud 1996). If possible, we wished to apply the same continuum model to both LMC X-4 and SMC X-1 spectra to allow for a direct comparison. We tested continuum models including Negative and Positive EXponential (NPEX, e.g., Mihara et al. 1998), a power law with a Fermi–Dirac cutoff (FDCut, Tanaka 1986), and a power law with a high-energy cutoff (White et al. 1983). We found that the FDCut and high-energy cutoff had slightly higher reduced χ² values and large residuals at high energies. For these reasons, we selected NPEX as our best continuum model. In Xspec, our NPEX model was defined as

$$\begin{eqnarray*}&&f(E)={n}_{1}({E}^{-{\alpha }_{1}}+{n}_{2}{E}^{-{\alpha }_{2}}){e}^{-E/{kT}},\end{eqnarray*}$$

where we fixed α₂ = −2.

In addition to NPEX, our spectral model also included an absorbing column (tbnew), a blackbody with kT ∼ 0.17 keV, and several Gaussian emission lines at 6.4 keV (Fe Kα), 1.02 keV (Ne x Lyα), 0.91 keV (Ne ix), and 0.65 keV (O viii Lyα). Each of these emission lines has been previously detected in LMC X-4 spectra with the Chandra High Energy Transmission Grating Spectrometer (Neilsen et al. 2009) and in SMC X-1 spectra with the Chandra ACIS instrument (Vrtilek et al. 2001, 2005). To reduce degeneracy in the blackbody model components, we fixed the widths of the Ne x Lyα, Ne ix, and O viii Lyα lines to the values found by Neilsen et al. (2009). In all observations, we found that the Fe Kα line was quite broad and that a 0.5 keV line width provided a good fit. However, in Observation S4 the spectrum also required a narrow (0.1 keV) component, as also seen by Neilsen et al. (2009).

To reduce degeneracies between the absorption and the blackbody component, we fixed the absorbing column density to the Galactic value in the direction of our sources, which we calculated using the HI4PI Map (HI4PI Collaboration et al. 2016) via the HEASARC ${N}_{{\rm{H}}}$ calculator. These values were 1 × 10²¹ cm⁻² for LMC X-4 and 3 × 10²¹ cm⁻² for SMC X-1.

For the LMC X-4 spectra, we used the elemental abundances described in Hanke et al. (2010) to account for the LMC's lower metallicity relative to Galactic abundances. The SMC X-1 spectral fits were performed using abundances from Wilms et al. (2000). For both sources, we used the cross sections from Verner et al. (1996).

The phase-averaged spectra and the residuals to the model fit for both data sets are shown in Figures 6 and 7. The spectral parameters and their uncertainties are given in Tables 5 and 6. We chose not to model XMM-Newton and NuSTAR in overlapping energy ranges so as to improve the model fit by minimizing differences in the response functions from these two observatories.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  LMC X-4 Phase-averaged Spectral Parameters^a

Parameter	Observation L1	Observation L2	Observation L3	Observation L4
${{kT}}_{\mathrm{BB}}$ (keV)	0.168 ± 0.006	0.168 ± 0.006	0.161 ± 0.004	0.160 ± 0.006
ABB (keV)	(5.2 ± 0.3) × 10−4	(7.0 ± 0.4) × 10−4	(2.27 ± 0.09) × 10−4	(3.2 ± 0.2) × 10−4
α1	0.55 ± 0.03	0.41 ± 0.02	0.55 ± 0.06	0.66 ± 0.04
${A}_{{\alpha }_{2}}$	(2.3 ± 0.1)$\,\times \,{10}^{-3}$	(4.2 ± 0.1) × 10−3	(5.5 ± 0.3) × 10−3	(2.9 ± 0.1) × 10−3
kTfold (keV)	6.1 ± 0.1	6.26 ± 0.05	5.81 ± 0.07	6.21 ± 0.07
log10(${F}_{3-40\mathrm{keV}}$)	−9.49 ± 0.01	−9.38 ± 0.01	−9.97 ± 0.02	−9.66 ± 0.01
E${}_{\mathrm{Fe}{\rm{K}}\alpha }$ (keV, fixed)	6.4	6.4	6.4	6.4
${\sigma }_{\mathrm{Fe}{\rm{K}}\alpha }$ (keV, fixed)	0.5	0.5	0.5	0.5
${A}_{\mathrm{Fe}{\rm{K}}\alpha }$	(3.3 ± 0.3) × 10−4	(2.2 ± 0.2) × 10−4	(1.04 ± 0.09) × 10−4	(1.3 ± 0.2) × 10−4
${E}_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$ (keV, fixed)	1.02	1.02	1.02	1.02
${\sigma }_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.003	0.003	0.003	0.003
${A}_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$	(7 ± 1) × 10−4	(4 ± 2) × 10−4	(1.2 ± 0.4) × 10−4	(1.9 ± 0.7) × 10−4
${E}_{\mathrm{Ne}{\rm{IX}}}$ (keV, fixed)	0.91	0.91	0.91	0.91
${\sigma }_{\mathrm{Ne}{\rm{IX}}}$ (keV, fixed)	0.003	0.003	0.22 ± 0.01	0.20 ± 0.01
${A}_{\mathrm{Ne}{\rm{IX}}}$	(3 ± 2) × 10−4	(1 ± 2) × 10−4	(8 ± 5) × 10−5	(2.2 ± 0.8) × 10−4
${E}_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.65	0.65	0.65	0.65
${\sigma }_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.003	0.003	0.003	0.003
${A}_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$	(1.0 ± 0.4) × 10−4	(1.1 ± 0.6) × 10−3	(5 ± 1) × 10−4	(7 ± 2) × 10−4
${c}_{\mathrm{EPIC} \mbox{-} \mathrm{pn}}$ (fixed)	1	1	1	1
cFPMA	2.10 ± 0.07	2.58 ± 0.08	3.0 ± 0.1	2.99 ± 0.09
cFPMB	2.13 ± 0.07	2.65 ± 0.09	3.1 ± 0.1	3.09 ± 0.09
${\chi }^{2}$	483.05	941.03	478.50	612.03
Degrees of freedom	440	908	425	568

Note.

^aFor the continuum model constant * tbnew * (cflux * npex + bbody + Gauss + Gauss + Gauss + Gauss). All errors are 90% confidence intervals.

Download table as:  ASCIITypeset image

Table 6.  SMC X-1 Phase-averaged Spectral Parameters^a

Parameter	Observation S1	Observation S2	Observation S4
kTBB (keV)	0.182 ± 0.001	0.179 ± 0.002	0.184 ± 0.001
ABB (keV)	(1.45 ± 0.01) × 10−3	(8.50 ± 0.09) × 10−4	(2.06 ± 0.01) × 10−3
${\alpha }_{1}$	0.402 ± 0.007	0.44 ± 0.02	0.397 ± 0.007
${A}_{{\alpha }_{2}}$	(1.01 ± 0.07) × 10−3	(1.9 ± 0.2) × 10−3	(1.18 ± 0.08) × 10−3
${{kT}}_{\mathrm{fold}}$ (keV)	5.64 ± 0.07	5.2 ± 0.1	5.56 ± 0.07
log10(${F}_{3-40\mathrm{keV}}$)	−9.260 ± 0.003	−9.450 ± 0.006	−9.119 ± 0.003
E${}_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{broad}}$ (keV, fixed)	6.4	6.4	6.4
${\sigma }_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{broad}}$ (keV, fixed)	0.5	0.5	0.5
${A}_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{broad}}$	(3.3 ± 0.5) × 10−4	(2.4 ± 0.3) × 10−4	(3.6 ± 0.6) × 10−4
E${}_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{narrow}}$ (keV, fixed)	N/A	N/A	6.4
${\sigma }_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{narrow}}$ (keV, fixed)	N/A	N/A	0.1 (fixed)
${A}_{\mathrm{Fe}{\rm{K}}\alpha ,\mathrm{narrow}}$	N/A	N/A	(2 ± 3) × 10−5
${E}_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$ (keV, fixed)	1.02	1.02	1.02
${\sigma }_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.003	0.003	0.003
${A}_{\mathrm{Ne}{\rm{X}}\mathrm{Ly}\alpha }$	(2.2 ± 0.7) × 10−4	(1.3 ± 0.5) × 10−4	(2.0 ± 0.8) × 10−4
${E}_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.65	0.65	0.65
${\sigma }_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$ (keV, fixed)	0.003	0.003	0.003
${A}_{{\rm{O}}{\rm{VIII}}\mathrm{Ly}\alpha }$	(3.1 ± 0.2) × 10−3	(1.7 ± 0.5) × 10−3	(4.3 ± 0.3) × 10−3
${c}_{\mathrm{EPIC} \mbox{-} \mathrm{pn}}$ (fixed)	1	1	1
cFPMA	3.21 ± 0.02	3.40 ± 0.05	2.63 ± 0.02
cFPMB	3.26 ± 0.03	3.48 ± 0.05	2.68 ± 0.02
${\chi }^{2}$	1317.43	725.34	1452.61
Degrees of freedom	963	614	973

Note.

^aFor the continuum model constant * tbnew * (cflux * npex + bbody + Gauss + Gauss + Gauss + Gauss). All errors are 90% confidence intervals.

Download table as:  ASCIITypeset image

For all LMC X-4 and SMC X-1 spectra, we fit the models jointly to the XMM-Newton, NuSTAR FPMA and FPMB spectra. The parameters in the NuSTAR spectra are tied to those in the XMM-Newton spectrum via a cross-calibration constant that accounts for differences in observed flux between the telescopes. The constants for NuSTAR FPMA and FPMB are in good agreement in all spectra. However, as can be seen in Tables 5 and 6, there is a discrepancy between the XMM-Newton and NuSTAR calibration constants, with the NuSTAR constants being 2–3 times the XMM-Newton constant (when ${c}_{{XMM}}=1$). This issue was even more pronounced for the SMC X-1 spectra, where XMM-Newton was in timing mode. We investigated the cross normalization in detail and found that our choice of wide XMM-Newton source extraction regions and non-overlapping energy when fitting drove the cross normalization unrealistically high. We verified that changes in source extraction region did not impact the spectral or pulse profile shapes. We modeled the joint NuSTAR and XMM-Newton spectra in the overlapping 3–10 keV range and found that the NuSTAR and XMM-Newton fluxes in this energy range agreed within 10%. A full exploration of the cross normalization is outside the scope of this work; however, we are confident that the values shown in Tables 5 and 6 are a reflection of our analysis steps and do not reflect the relative fluxes measured by XMM-Newton and NuSTAR.

2.3.2. Phase-resolved Spectroscopy

The pulse profiles of LMC X-4 and SMC X-1 shown in Figures 4 and 5 show changes in shape and phase over the course of a single superorbital cycle.

For LMC X-4, the hard (8–60 keV) and soft (0.5–1 keV) pulse profiles from Observations L1 and L2 are out of phase. In Observation L1 the soft pulses are slightly less than 180° out of phase, while in Observation L2 they appear to be closer to 180° out of phase. By contrast, the hard and soft pulse profiles in Observation L3 are almost completely in phase. The pulse profiles, and in particular the hard pulses, in Observation L4 are weakly detected due to the pulse dropout phenomenon that occurred during this observation (see Brumback et al. 2018). Despite this, we observe that the hard and soft pulsations appear out of phase.

Independent of pulse dropout behaviors, we also observe changes in pulse shape with superorbital phase in LMC X-4. The hard pulse profiles in Observations L1 and L3 are relatively smooth single peaks, while the hard pulse profile in Observation L2 has become broad and flat. In general, the soft pulse shapes in all LMC X-4 observations are rounded single peaks. We also observe a shift in relative strength between hard and soft pulsations in Observation L2, apparently driven by a change in the hard pulsed fraction.

The pulse profile for SMC X-1 is double-peaked. With the energy-resolved pulse profiles for SMC X-1, we find that the profiles for Observations S1 and S4 are extremely consistent in shape and relative phase; both hard and soft profiles are in phase with each other and both show two peaks of approximately equal strength in both the hard and soft pulses. These shapes are different than those seen in Observation S3, where the hard pulses show one strong and one weak peak, while the soft pulses have merged into a broad single peak.

Because the hard pulsations are caused by the pulsar beam and the soft pulsations originate from accretion disk reprocessing (e.g., Hickox et al. 2004), the consistency in pulse phase between observations from the same superorbital phase, particularly Observations S1 and S4, shows that we have observed a complete precession cycle of the inner accretion disk.

Our spectroscopic analysis of these data indicates that the broadband X-ray spectra of LMC X-4 and SMC X-1 are well described by an absorbed power law and a soft blackbody component. In both sources, the blackbody temperature changes very little with superorbital phase. In SMC X-1, we also find very little variation in the α₁ parameter with superorbital phase. There is, however, some variation in the strength of the second power law (${A}_{{\alpha }_{2}}$), which indicates that the overall shape of the hard continuum is changing slightly with superorbital phase. In LMC X-4, we observe changes in power-law shape through variation in both α₁ and ${A}_{{\alpha }_{2}}$.

We would expect that Observations S1 and S4 would have generally the same spectral shape since they were taken at the same superorbital phase, and the same applies to Observations L1 and L4. We do find good agreement between the spectral parameters in Observations S1 and S4, where the only notable differences are a slightly stronger second power-law normalization and the presence of a narrow Fe Kα feature in Observation S4. We find less good agreement between the spectral parameters in Observations L1 and L4; however, we do not consider these discrepancies to be problematic in view of the fact that these two observations sample different accretion and pulse behaviors (Brumback et al. 2018). In Observation L1, our phase-averaged spectrum reflects a strongly pulsed time interval between bright accretion flares, whereas in Observation L4 our spectrum is drawn from a weakly pulsed pre-flare interval. The hardness ratios vary between these two states, implying that the shape of the spectrum changes (see Figure 1 in Brumback et al. 2018). In Brumback et al. (2018) we suggest that these different pulse behaviors could be driven by changing emission geometries during the accretion flares. If this is indeed the case, we would expect to see differences in the spectral shape during this process.

In our phase-resolved analysis, the parameters allowed to vary within each spectrum were the blackbody normalization, the overall flux of the power law, the primary power-law index α₁, the secondary power-law normalization, the power-law folding energy, and the Fe Kα line normalization. Across the phase-resolved spectra for both LMC X-4 and SMC X-1, we only find clear, coherent changes with pulse phase in the blackbody normalization, the power-law flux, and the Fe Kα line normalization. The other parameters (α₁, folding temperature, and second power-law normalization) either are consistent with being constant or show variations that are difficult to describe physically due to degeneracies within the model and reduced signal-to-noise ratio in the spectra.

An example of the smooth variations in power-law flux and blackbody normalization is shown in Figure 8 for Observation L3, where the variation in power-law flux is overplotted with the hard NuSTAR pulse profile, and the blackbody normalization is overplotted with the soft XMM-Newton pulse profile. The other LMC X-4 and SMC X-1 observations show similarly good agreement between these two spectral parameters and the pulse profiles, and so for the sake of brevity we do not show them here. The agreement of these parameters with their respective pulse profiles is significant because it indicates that the pulse profiles (measured in count rates) are a reasonable proxy for the spin-resolved power law and blackbody flux. This agreement allows us to directly fit the energy-resolved NuSTAR and XMM-Newton pulse profiles in our warped disk model and to assume that the pulse profiles are following the changes in strength of the power law and blackbody.

Zoom In Zoom Out Reset image size

Hickox et al. (2004) found that disk reprocessing is a ubiquitous feature of bright X-ray pulsars. HV05 used a simple warped disk model to describe the differences in shape and phase between the hard and soft pulsations in SMC X-1 as they vary across the superorbital cycle. Hung et al. (2010) used the same model to qualitatively describe the pulse profiles in LMC X-4. While these previous works demonstrated the success of the HV05 disk model, neither used observations within a single disk precession cycle to examine the periodicity of the disk and determine whether this model can describe pulse behavior over a complete disk cycle.

We seek to verify whether the HV05 model can reproduce the changes in pulse shape and phase seen over a complete disk precession cycle in both LMC X-4 and SMC X-1. The warped disk model used in this analysis is the same as that presented by HV05 and used by Hung et al. (2010).

HV05 describes the warped inner region of the accretion disk as a series of concentric circles that are inclined and rotated relative to each other. This geometry is based on the well-constrained disk of the bright X-ray binary Her X-1 (Scott et al. 2000; Leahy 2002). The precise geometry of this disk is set by the radii of the inner and outer circles and their respective inclination angles (r_in, r_out, θ_in, θ_out). We provide a schematic diagram of the HV05 disk and beam geometries in Figure 9.

Zoom In Zoom Out Reset image size

We also use two simple beam geometries, a pencil and a fan beam, to model the neutron star's beam geometry. In both cases, the beams are modeled as two-dimensional Gaussians with width σ_b. The location of the beam on the surface of the neutron star is defined by the angle out of the plane of rotation and the azimuthal angle (θ_b, ϕ_b). We define the coordinate system such that the poles align with the rotation axis and θ = 0 lies along the equator, and the neutron star's rotation is parallel to the disk axis. In the fan beam model, we also define a beam opening angle (θ_fan), which is set to 0 in the pencil model. For simplicity, the fan beam model is a pure fan beam, without an embedded pencil beam.

In this model, the observer is set at a fixed angle (${\theta }_{\mathrm{obs}}$), which, if the neutron star rotates within its orbital plane, is related to the inclination angle for the system by i = 90° − θ_obs. The beam pattern is then rotated and the regions of the disk visible from the neutron star are illuminated. The disk is assumed to be opaque and it immediately reradiates the absorbed emission as a blackbody spectrum. This assumption requires that the light-crossing time and disk cooling time be shorter than the pulse period of the neutron star. The light-crossing time for a disk surface of approximately 10⁸ cm is ∼10 ms. For the cooling time, Endo et al. (2000) suggested that this timescale can be estimated as the thermal energy of the disk divided by the luminosity. For general parameters such as a Compton thick disk, a blackbody temperature of ${{kT}}_{\mathrm{BB}}=0.18\,\mathrm{keV}$, and a soft X-ray luminosity of 10³⁷ erg s⁻¹, HV05 estimate the cooling time as ∼10⁻⁵ s. Both the light-crossing time and the cooling time are shorter than the pulse periods of LMC X-4 and SMC X-1, and therefore we assume immediate reprocessing by the disk.

Emission seen by the observer is calculated at 30 pulse phases and eight equally spaced disk phase intervals, where disk phase zero is defined as when the neutron star first emerges from behind the disk, consistent with the start of the superorbital high state. For each beam geometry and disk rotation phase the luminosity of the beam and the luminosity of the disk regions visible to the observer are calculated, and simulated hard (beam) and soft (disk) pulse profiles are made. The HV05 model does not include the effects of light bending on the emission viewed by the observer.

In our model we constrain the disk surface between an inner radius of 0.8 × 10⁸ cm and an outer radius of 1 × 10⁸ cm, which HV05 found reproduced the observed SMC X-1 blackbody temperature. We initially set the observer angle to 20° because this agrees with estimates of the orbital inclination of ∼70° for both SMC X-1 and LMC X-4 (Reynolds et al. 1993; van der Meer et al. 2007). While this value worked well for the SMC X-1 models, we found that we could not reproduce the pulse profiles of Observation L3 with an observer angle of 20°. We tested a range of observer angles from 5° to 40° and found that the pulse profiles of Observation L3 could only be reproduced with an observer angle of 40°. We set the outer disk angle to be 45° for both sources; this is within the range of disk inclination of 25°–58° estimated for SMC X-1 by Lutovinov et al. (2004), and we found this angle necessary to reproduce the observed pulse profiles of LMC X-4. Our outer disk angle also agrees with hydrodynamic simulations that Larwood et al. (1996) used to find stable precession in tilted accretion disks with outer disk angles of 45°. We fixed the inner disk angle to a smaller value of 10°. We found that a beam half-width of 30° fit all observations well. While the disk geometry was allowed to vary between LMC X-4 and SMC X-1, we used the same disk parameters to describe the observations from each source, but allowed the beam parameters to change between observations.

When fitting the pulse profiles to data, we allow the overall intensity of the simulated pulses to vary so that the intensity matches that of the observed hard pulsations.

To simulate pulse profiles for both the pencil and fan beam models of LMC X-4, we began fitting pulse profiles with the brightest observation in the data set: Observation L3. We first simulated the hard pulse profile shape to match the observed data and then adjusted the disk parameters until we found reasonable agreement in the soft pulse profiles. We then kept the disk parameters the same for the other three LMC X-4 observations and varied the beam height and azimuth (${\theta }_{{\rm{b}}}$, ϕ_b), which was necessary to match the other pulse shapes in the observation series. We note that these changes in beam location do not necessarily represent physical changes in the accretion column, but rather reflect the varying effects of light bending or other phenomena not included in the HV05 model. For each observation, we allowed the disk to precess and calculated pulse profiles for each precession phase. We fit the three SMC X-1 observations in the same way. The best-fit parameters for both the pencil and fan beam configurations are listed in Table 7.

Table 7.  Disk Model Parameters

	LMC X-4		SMC X-1
Parameter	Pencil Beam	Fan Beam	Pencil Beam	Fan Beam
rin (108 cm)	0.8	0.8	0.8	0.8
rout (108 cm)	1	1	1	1
Inner tilt θin (deg)	10	10	10	10
Outer tilt θout (deg)	45	45	45	45
Twist angle ϕtw (deg)	−130	−130	−130	−130
Beam1 angle from rotational plane θb1 (deg)	60a,b,c, 75d	60b, 70a,c, 75d	60e,f, −50g	60e,f, −40g
Beam2 angle from rotational plane ${\theta }_{{\rm{b}}2}$ (deg)	60a,c, 65d, −60b	60d, 70a,c, −60b,	60e,f,g	60e,f,g
Beam1 azimuth ϕb1 (deg)	0	0	0	0
Beam2 azimuth ϕb2 (deg)	110d, 130a,c, 160b	110d,120a,c, 160b	180e,f, 185g	180e,f,g
Beam half-width σb (deg)	30d, 45a,c, 60b	30	60	30
Fan beam opening angle θfan (deg)	0	15d, 20b, 25a,c	0	30e,f,g
Observer elevation θobs (deg)	40	40	20	20

Notes.

^aThis value required for Observation L2. ^bThis value required for Observation L3. ^cThis value required for Observation L4. ^dThis value required for Observation L1. ^eThis value required for Observation S1. ^fThis value required for Observation S4. ^gThis value required for Observation S2.

Download table as:  ASCIITypeset image

To find the best fit to the soft pulses, we estimated the goodness of fit between the simulated pulse profiles produced at different precession phases and the observed pulse profile by calculating $r=\sum ({P}_{\mathrm{obs}}({\phi }_{\mathrm{spin}})-{P}_{\mathrm{sim}}({\phi }_{\mathrm{spin}}))/\overline{{P}_{\mathrm{obs}}}$, where P_obs is the observed pulse profile and P_sim is the simulated pulse profile, and identifying the disk phases with the lowest r value. These best-fit disk phases are highlighted in green in Figures 10 and 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In LMC X-4, the HV05 model is able to describe the shape of the hard pulsations with the exception of Observation L4, which has extremely weak pulsations. The lack of significant pulsations in this observation is possibly due to pulsation dropout in association with super-Eddington accretion flares, and the timing properties of this observation are discussed in Brumback et al. (2018). For the purposes of this analysis, the effect of weak pulsations in Observation L4 results in poor constraints on the beam profile.

By allowing the disk to precess, the HV05 model successfully reproduces the shape of most of the LMC X-4 soft pulsations in at least one disk phase. However, the HV05 model struggles to reproduce the phase of the soft pulsations in Observation L3 (which are nearly in phase with the observed hard pulsations). This is most likely because of the broad beam parameters necessary to create single-peaked pulse profiles. We indicate the best-fit disk phases for the fan beam configuration in Figure 10, but we note that this is likely not a valid constraint on the disk precession phase. The pencil beam configuration produced similar results to those shown in Figure 10, and so we do not include these figures for the sake of conciseness.

We also find good fits to the observed hard pulse profiles for SMC X-1. We show the results of the fan beam configuration in Figure 11, and again do not show the similar results from the pencil beam configuration for the sake of space. We found that the HV05 model struggled to reproduce the soft pulses observed in Observation S2. The challenges in simulating these pulse profiles likely arise from the hard pulsations having a double-peaked profile; when the hard profile was double-peaked the model strongly preferred a soft profile that was double-peaked as well. We found that for a double-peaked hard profile, the HV05 model was not able to return a single-peaked soft profile as broad as the observed profile.

Despite challenges to modeling presented by Observations L3 and S2, we find that in both LMC X-4 and SMC X-1 the disk phase values corresponding to our best-fit soft pulse profiles are consistent with a complete precession cycle of the inner accretion disk (Figures 12 and 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our simulation of the hard and soft profiles confirmed general conclusions drawn by Hickox & Vrtilek (2005), including that the pulse profile shape is more dependent on the beam geometry than the disk geometry, and that the double- and single-peaked pulse profiles seen in these sources strongly prefer non-antipodal beam geometry. This preference can be seen in Table 7, where negative values of θ_b1 and θ_b2 are only found in Observations L3 and S2.