The Stochastic X-Ray Variability of the Accreting Millisecond Pulsar MAXI J0911–655

The X-ray transient MAXI J0911 was discovered on 2016 February 19 (Serino et al. 2016) with the MAXI/GSC. The source was immediately associated with the globular cluster NGC 2808, a position that was later confirmed by Swift/XRT (Kennea et al. 2016) and Chandra (Homan et al. 2016) observations. Subsequent monitoring with INTEGRAL and Swift/XRT has shown that MAXI J0911 has remained active, showing a persistent flux of about 7 mCrab (Ducci et al. 2016). At the time of writing, this source had yet to transition into quiescence, placing the duration of the outburst at approximately 1 yr.

The nature of the compact object in MAXI J0911 was settled when its 340 Hz pulsations were discovered by Sanna et al. (2017). Using pointed XMM-Newton and NuSTAR observations, these authors studied the pulsations and showed that the AMXP is set in a compact binary with a 44.3 minute orbital period and a companion star of $\gt 0.024$ ${M}_{\odot }$. The nature of the companion is not definitively constrained, but it is likely either a hot, helium white dwarf or an old brown dwarf. Additionally, Sanna et al. (2017) found that the energy spectrum of MAXI J0911 is relatively hard, with a ${\rm{\Gamma }}=1.7$ power law dominating over a weak ${kT}\sim 0.5$ blackbody component.

In this work, I analyze the epic-pn data of two XMM-Newton observations of MAXI J0911. The first observation took place on 2016 April 24 (ObsID 0790181401) and the second on 2016 May 22 (ObsID 0790181501). For both observations, the epic-pn camera was operated in timing mode, yielding event list data at a time resolution of 29.56 $\mu {\rm{s}}$.

To obtain science-grade products, the data were processed with sas version 15.0.0 using the most recent calibration files available. Standard data-screening criteria were applied to the data, selecting only those events in the $0.4\mbox{--}10\,\mathrm{keV}$ energy range with pattern ≤ 4 and screening flag = 0. The source data were extracted from a 15-bin-wide rectangular region with the rawx coordinates [31:45]. Finally, the sas tool barycen was applied to correct the event arrival times to the solar system barycenter using the source coordinates of Homan et al. (2016).

The background estimate was obtained similarly but from a 3-bin-wide region with the rawx coordinates [3:5]. Since this region is smaller than the source extraction region, I multiplied the resulting background count rate by a factor of 5 to ensure that both source and background rates reflected a comparable effective area.

The source count rates obtained are 40(1) and 35(2) counts s⁻¹ for the first and second observations, respectively. The associated background count rates are 0.8(4) and 0.4(2) counts s⁻¹, respectively.

NuSTAR observed MAXI J0911 on 2016 May 24 (ObsID 90201024002) and again 180 days later on 2016 November 23 (ObsID 90201042002). The data were processed using the nustardas software pipeline version 1.6.0. This was done separately for each of the two focal plane modules (FPMA and FPMB).

The source events in the $3\mbox{--}79\,\mathrm{keV}$ energy band were extracted from a circular region with a 50'' radius that was centered on the image source position. The filtered event arrival times were then corrected to the solar system barycenter using the barycorr tool, again based on the source coordinates of Homan et al. (2016). The background events were extracted similarly but with a source region centered in the background field. The resulting source count rates are 2.0(2) and 3.0(2) counts s⁻¹ for the first and second observations, respectively, with a background count rate of 0.01(1) counts s⁻¹ in both cases.

A stochastic timing analysis is performed on the cleaned and processed X-ray data. First, all the arrival times are adjusted for the binary motion of the neutron star based on the ephemeris of Sanna et al. (2017). The event lists are then binned into ∼250 s light curves using a time resolution of 480 $\mu {\rm{s}}$ (the length of the time series is adjusted so that the number of time bins is a power of 2). For each of these light curves, I compute the Fourier transform, from which a Leahy-normalized power spectrum is calculated (Leahy et al. 1983). These power spectra have a frequency resolution of ∼0.004 Hz and a limiting Nyquist frequency of 1040 Hz. All individual power spectral estimates are averaged per ObsID and corrected for Poisson noise.

For the XMM-Newton data, the Poisson noise spectrum is modeled as a constant. That noise power is measured between 400 and 1000 Hz, where no signal features are observed, and then subtracted from the averaged power density spectrum.

In the case of the NuSTAR data, the instrument dead time needs to be taken into account, as it has a comparatively long timescale of ∼2.5 ms (Harrison et al. 2013). An elegant way of dealing with the NuSTAR dead time was proposed by Bachetti et al. (2015) and involves computing the complex-valued cross spectrum from the two FPM light curves. The real component of the cross spectrum, also known as the co-spectrum, can then be adopted as an estimator of the source power density spectrum. Because, to good approximation¹ , the dead time is independent between detectors while the signal itself is wholly correlated, this approach effectively filters out any dead-time modulation without the need to model it.

I would point out, however, that the gain achieved through considering the co-spectrum does not come free. The cost of correlating the two FPM light curves is that the resulting signal-to-noise ratio is lower than what would have been obtained if the two light curves were simply added in the time domain. Specifically, the uncertainties on the co-spectral powers are larger than those obtained from a regular power spectrum by up to a factor of $\sqrt{2}$ (see Appendix A for details).

Given that the observed event rate (0.5 s count⁻¹) is much larger than the timescale of the dead time, the effect of the dead time will be relatively small, and the power spectrum will only be sensitive at comparatively low frequencies. I would therefore argue that, in this case, it is more appropriate to add the two FPM light curves for the increased sensitivity, rather than correlating them to reduce the dead-time modulation.

The NuSTAR power spectra are constructed by adding the concurrent FPM data. For each light curve, the Fourier transforms are computed and normalized and subsequently averaged in the complex domain. From this averaged Fourier transform, a Leahy-normalized power spectrum is constructed. Note that this power spectrum will have an expected Poisson noise level of 1 (see Appendix A). All segments of the ObsID are averaged to a single power spectrum. This power spectrum shows some dead-time modulation, which is modeled heuristically by taking all powers above 25 Hz and fitting the scaled sinc function

$$\begin{eqnarray*}&&f(\nu | A,B,{t}_{D})=A-2{{Bt}}_{D}{\rm{sinc}}(2\pi \nu {t}_{D}),\end{eqnarray*}$$

where A and B are a mean and a scale parameter, respectively, and t_D is the characteristic dead time. For the first NuSTAR observation, the best fit gives ${\chi }^{2}/\mathrm{dof}=147/158$ with parameters $A=1\pm 1\times {10}^{-4}$, $B=1.5\pm 0.1$, and ${t}_{D}=2.56\pm 0.06$ ms (Figure 1). Below 10 Hz, this curve is near constant at a level of ${P}_{\mathrm{noise}}\simeq A-2{{Bt}}_{D}=0.992$, which I adopt as the Poisson noise level and subtract from the data. The second NuSTAR observation has similar results and differs only in the scale, $B=2.32\pm 0.14$, so that the Poisson noise level is marginally higher, ${P}_{\mathrm{noise}}\simeq 0.994$.

Zoom In Zoom Out Reset image size

Finally, the Poisson-noise–corrected averaged power spectra are renormalized to squared fractional rms and subsequently described in terms of a multi-Lorentzian model (Belloni et al. 2002). Each component, $L(\nu | r,Q,{\nu }_{\max })$, is a function of Fourier frequency ν and characterized by three parameters. Here ${\nu }_{\max }={\nu }_{0}\sqrt{1+1/4{Q}^{2}}$ is the characteristic frequency, and $Q={\nu }_{0}/W$ is the quality factor, where ${\nu }_{0}$ gives the centroid frequency and W is the full width at half maximum. The amplitude of an individual component is expressed in terms of fractional rms, r, defined as

$$\begin{eqnarray*}&&{r}^{2}=P={\int }_{0}^{\infty }L(\nu )d\nu ,\end{eqnarray*}$$

where P is the integrated power. A component is considered to be significant if its integrated power has a single trial significance greater than 3; that is, if $P/{\sigma }_{P}\geqslant 3$.

The overall morphology of the power spectrum is reminiscent of an extreme island state atoll source (van der Klis 2006). For such a state, the three noise terms may be identified as the break, hump, and ${\ell }$ow components, each having a progressively higher frequency. The two QPOs, seen atop the hump component, may then be identified as the LF QPO and its harmonic (LF₂; van Straaten et al. 2003; Altamirano et al. 2005). The energy spectrum of MAXI J0911 is relatively hard (Sanna et al. 2017), which is consistent with an atoll source in this source state. The power spectral components, as shown in Table 1, have therefore been labeled according to this terminology.

A notable difference between MAXI J0911 and the wider class of atoll sources is that the frequencies measured in this work are about an order of magnitude lower than what is usually observed (see, e.g., van Doesburgh & van der Klis 2017). Meanwhile, the observed fractional variability is remarkably high.

While extreme, these properties are in line with the source-state evolution seen in atoll sources. As the energy spectrum increases in hardness, the power spectral features are observed to shift to lower frequencies while broadening and increasing in amplitude (van der Klis 2006). In this context, the source state of MAXI J0911 would represent an outlier along the evolutionary track, suggesting that this system is accreting in a geometry not normally accessible by regular atoll-type objects.

The identification of power spectral features in neutron star sources is usually confirmed by considering the frequency relations as a function of the highest-frequency component in the spectrum. However, with the high-frequency end of the power spectrum being poorly constrained, this approach is not possible for MAXI J0911. Instead, I test the frequency relation of the observed features with respect to the hump component. In atoll sources, these features have nearly fixed frequency relations (van Straaten et al. 2002; Altamirano et al. 2008a), so that their frequency ratios are approximately constant (van Doesburgh & van der Klis 2017). Hence, the frequency ratios of the measured components relative to the hump component frequency are shown in Figure 4. While the absolute frequencies in MAXI J0911 are low, the frequency ratios are consistent with those in other atoll sources and AMXPs. This strongly supports the identification of these features as LF QPOs.

Zoom In Zoom Out Reset image size

An alternative interpretation of the QPOs may be that they are related to the mHz QPOs (Revnivtsev et al. 2001), which are attributed to marginally stable burning of accreted material on the neutron star surface (Heger et al. 2007). While the luminosity of MAXI J0911 is in the regime where such QPOs are observed, there are several issues with this interpretation. Both the QPO frequencies and the amplitudes in MAXI J0911 are higher by a factor of a few than is usual for mHz QPOs (Altamirano et al. 2008b; Linares et al. 2012), and so far no type I X-ray burst has been observed in MAXI J0911. More importantly, the QPO amplitudes show a steep increase as a function of energy, strongly disfavoring a scenario in which they originate from the neutron star surface.

Yet another possibility might be that periodic obscuration along the line of sight is giving rise to these QPOs. However, again, this seems unlikely. Dipping QPOs tend to have higher frequencies and amplitudes than what is observed in MAXI J0911 (Homan et al. 1999). Furthermore, the rms amplitude of a dipping QPO is expected to be independent of energy, which does not correspond with the observations.

There is only one other neutron star system that shows a power spectrum that is similar to the one observed in MAXI J0911. That source is the 599 Hz accreting millisecond pulsar IGR J00291+5934 (IGR J00291; Linares et al. 2007).

Observed with RXTE in the $2\mbox{--}20\,\mathrm{keV}$ band, the power spectrum of IGR J00291 shows a $0.01\mbox{--}100\,\mathrm{Hz}$ integrated power of $40 \% \mbox{--}60 \%$ rms (Linares et al. 2007). The break, hump, and low components all have comparable frequencies and amplitudes to those seen in MAXI J0911. Additionally, two very-low-frequency QPOs have been detected in IGR J00291. With frequencies of 22 and 44 mHz, however, they sit atop the break component rather than the hump component and thus cannot be identified as LF QPOs like those in MAXI J0911.

A very prominent 8 mHz flaring has also been detected in IGR J00291 (Ferrigno et al. 2017). This feature is not seen MAXI J0911, but if it is due to "heartbeat" variability (Belloni et al. 2000; Altamirano et al. 2011), as suggested by Ferrigno et al. (2017), then such flaring may be a transient phenomenon, and MAXI J0911 may in fact be a good candidate to monitor for such variability.

The similarity of their broadband power spectra can be taken to indicate that both MAXI J0911 and IGR J00291 are accreting in a comparable source state. That source state is then characterized by a highly energetic Comptonizing medium and a large fractional variability. As noted by Linares et al. (2007), these properties are naturally explained if the optically thick accretion disk is truncated at a relatively large radius. An extended inner flow would then set the high fractional amplitude and the energy dependence, while the comparatively cool outer disk would set the low dynamical frequency (Churazov et al. 2001; Gilfanov et al. 2003). How such a configuration would interact with the stellar magnetic field remains unclear (see Patruno et al. 2016 for a discussion).

There is no obvious observable parameter indicating why these particular sources would have a larger truncation radius or, more generally, such an extreme source state. Both systems are millisecond pulsars, but they have somewhat different properties (Patruno 2016; Sanna et al. 2017). Their spin frequencies differ by nearly a factor of 2. With an orbital period of 2660 s, MAXI J0911 is in a tighter orbit than IGR J00291, which has a period of 8844 s. The mass limits on their companions are similar, and both have accretion luminosities of ∼10³⁶ erg s⁻¹ (Falanga et al. 2005; Sanna et al. 2017). However, the same can be said for other accreting millisecond pulsars. For instance, SAX J1808.4–3658 is more similar to IGR J00291 than MAXI J0911, and yet the power spectra of SAX J1808.4–3658 are far more consistent with atoll sources at those luminosities (Bult & van der Klis 2015).

Another peculiar property of MAXI J0911 is the duration of its outburst. Unlike IGR J00291, which shows a few week-long outbursts every 4–5 yr, MAXI J0911 has been continuously accreting for at least 1 yr. From the two NuSTAR observations it follows that, during this time, the source state has remained unchanged. The only other AMXP to have shown such a long outburst is the 377 Hz intermittent pulsar HETE J1900.1–2455 (Kaaret et al. 2006), which remained active for about 10 yr. The power spectrum of that source, however, is that of a regular (extreme) island state atoll source (Papitto et al. 2013). This would suggest that neither intermittency nor outburst duty cycle has much bearing on the stochastic variability of these pulsars.

Potentially, the accretion geometry of MAXI J0911 and IGR J00291 can be accounted for by the detailed properties of the neutron star, depending on the magnetic field strength, stellar mass, and spin and magnetic alignment angles. While such properties are difficult to measure, detailed pulse profile modeling for IGR J00291 and MAXI J0911 may provide further clues.

The LF QPOs observed in MAXI J0911 have the lowest frequency seen in any accreting neutron star system. Their frequency range is more similar to that of the LF QPOs observed in low-hard state black hole binaries (Klein-Wolt & van der Klis 2008). This is illustrated in Figure 4, which also shows frequency ratios for the black hole binary GX 339−4.

For black hole binaries, there is strong evidence that LF QPOs (or type C QPOs in black hole terminology) are caused by Lense–Thirring precession (Ingram et al. 2016, 2017). In this model, the accretion disk consists of two components: an inner hot accretion flow emitting a hard Comptonized spectrum and an outer thin disk emitting a softer thermal spectrum. If the black hole spin is misaligned with the accretion disk, frame-dragging effects apply a torque on the hot flow, causing it to precess as a solid body (Ingram et al. 2009). It is argued that this precession causes the LF QPOs. Qualitatively, this model is in good agreement with the power spectral features of MAXI J0911 and the hard energy spectra of their amplitudes.

If neutron star systems have a precessing inner accretion flow like that predicted for black holes, then their LF QPO frequencies are not straightforwardly related to the Lense–Thirring precession frequency (Altamirano et al. 2012; Bult & van der Klis 2015; van Doesburgh & van der Klis 2017). For instance, in neutron star systems, the stellar oblateness introduces an additional torque on the disk, giving rise to a retrograde precession term (Morsink & Stella 1999). For accreting millisecond pulsars, this situation is complicated further by the stellar magnetosphere, which truncates the accretion disk and has been proposed to introduce a third, magnetic torque (Lai 1999; Shirakawa & Lai 2002), leading to a second retrograde precession term.

Determining how, exactly, these various torques set the LF QPO frequency has the potential to yield constraints on how the magnetosphere interacts with the disk and on the neutron star properties themselves. The black hole–like power spectrum of MAXI J0911 makes this source an ideal target for such a study, possibly with the Neutron Star Interior Composition Explorer (Gendreau et al. 2012), which is scheduled for launch in 2017.

Consider a discrete time series x_j containing a total of N_γ counts. It is well known that, for a time series containing only Poisson noise, the elements of the Fourier transform, X_k, will be distributed as complex normal variates with zero mean and variance ${N}_{\gamma }/2$ (Leahy et al. 1983). If I let

$$\begin{eqnarray}&&{Z}_{k}=\sqrt{\displaystyle \frac{2}{{N}_{\gamma }}}{X}_{k}={A}_{k}{e}^{i{\phi }_{k}}+{Z}_{\mathrm{noise}}\end{eqnarray} \tag{ 1 }$$

be a normalized Fourier transform of a signal in the presence of Poisson noise, then A_k² gives the signal Leahy power, ${\phi }_{k}$ is its underlying phase, and ${Z}_{\mathrm{noise}}$ is a complex-valued noise term with independent standard normal variates in both its real and its imaginary parts.

Given an ensemble of M concurrent measurements, the Fourier transforms can be averaged coherently (dropping the frequency index k for convenience),

$$\begin{eqnarray}&&{\langle Z\rangle }_{M}={{Ae}}^{i\phi }+{\langle {Z}_{\mathrm{noise}}\rangle }_{M},\end{eqnarray} \tag{ 2 }$$

so that the variance reduces to $1/M$. By definition, the sum of squared standard normal variates follows a χ² distribution. Similarly, the sum of squared normal variates that have a nonzero mean, ${\mu }_{n}$, and unit variance will follow a noncentral χ² distribution, ${\overline{\chi }}^{2}(\nu ,\lambda )$, characterized by ν degrees of freedom and noncentrality,

$$\begin{eqnarray}&&\lambda =\displaystyle \sum _{n=1}^{\nu }{\mu }_{n}^{2}.\end{eqnarray} \tag{ 3 }$$

It then follows that the power of ${\langle Z\rangle }_{M}$ is distributed as ${\overline{\chi }}^{2}(2,{{MA}}^{2})$ scaled by a factor of $1/M$.

Given an ensemble of L coherently averaged powers, the combined averaged power is distributed as

$$\begin{eqnarray}&&{f}_{P}(y)=\displaystyle \frac{1}{{LM}}{\overline{\chi }}^{2}\left(\displaystyle \frac{y}{{LM}}| \nu =2L,\lambda ={{LMA}}^{2}\right),\end{eqnarray} \tag{ 4 }$$

which has a mean of

$$\begin{eqnarray}&&{\mu }_{P}=\displaystyle \frac{2}{M}+{A}^{2}\end{eqnarray} \tag{ 5 }$$

and a variance of

$$\begin{eqnarray}&&{\sigma }_{P}^{2}=\displaystyle \frac{4+4{{MA}}^{2}}{{{LM}}^{2}}.\end{eqnarray} \tag{ 6 }$$

Note that, for a signal that is only averaged incoherently (M = 1), Equation (4) reduces to the result derived by Groth (1975) for a signal in the presence of noise. For a pure noise signal (A = 0), these results reduce further to the familiar χ² distribution with $2L$ degrees of freedom and yield a Poisson noise level of 2 (van der Klis 1989).

For an ensemble of M concurrent time series (as in Equation (1)), there are $\tilde{M}=\left(\genfrac{}{}{0em}{}{M}{2}\right)=(M-1)M/2$ unique pairs, each of which can be correlated to give an estimate of the co-spectrum. For such a pair, say ${x}_{a}[j]$ and ${x}_{b}[j]$, the co-spectrum at a particular frequency index is given as

$$\begin{eqnarray}&&{C}_{{ab}}={\rm{Re}}[{Z}_{a}{Z}_{b}^{* }]={R}_{a}{R}_{b}+{I}_{a}{I}_{b},\end{eqnarray} \tag{ 7 }$$

where R_a and R_b are normal variates with mean $A\cos \phi$ and unit variance. The terms I_a and I_b are similarly distributed with a mean of $A\sin \phi$. It is useful to define the product ${Y}_{{ab}}={R}_{a}{R}_{b}$, for which the expectation and variance are given as

$$\begin{eqnarray}{\mathbb{E}}[{Y}_{{ab}}] & = & {\mathbb{E}}[{R}_{a}]{\mathbb{E}}[{R}_{b}]={A}^{2}{\cos }^{2}\phi ,\\ {\mathbb{V}}[{Y}_{{ab}}] & = & {\mathbb{E}}[{R}_{a}^{2}]{\mathbb{E}}[{R}_{b}^{2}]-{({\mathbb{E}}[{R}_{a}]{\mathbb{E}}[{R}_{b}])}^{2}\\ & = & 1+2{A}^{2}{\cos }^{2}\phi .\end{eqnarray} \tag{ 8 }$$

When averaging the product Y over all $\tilde{M}$ correlation pairs, the mean will reduce to this expectation. For the variance, it is important to realize that the pairs are not all independent. Specifically, the two products Y_ab and Y_ac have a covariance of

$$\begin{eqnarray}{\rm{Cov}}[{Y}_{{ab}},{Y}_{{ac}}] & = & {\mathbb{E}}[{R}_{a}^{2}{R}_{b}{R}_{b}]-{\mathbb{E}}[{Y}_{{ab}}]{\mathbb{E}}[{Y}_{{ac}}]\\ & = & {A}^{2}{\cos }^{2}\phi .\end{eqnarray} \tag{ 9 }$$

If I let ${{\rm{\Gamma }}}_{{nm}}$ be the full covariance matrix of the $\tilde{M}$ correlation pairs, then the variance of the averaged product term will be

$$\begin{eqnarray*}{\mathbb{V}}[{\langle Y\rangle }_{\tilde{M}}] & = & \displaystyle \frac{1}{{\tilde{M}}^{2}}\displaystyle \sum _{m=0}^{\tilde{M}-1}\,\displaystyle \sum _{n=0}^{\tilde{M}-1}{{\rm{\Gamma }}}_{{nm}}\\ & = & \displaystyle \frac{1}{{\tilde{M}}^{2}}\left\{\displaystyle \sum _{m=0}^{\tilde{M}-1}{{\rm{\Gamma }}}_{{mm}}+\displaystyle \sum _{m=0}^{\tilde{M}-1}\,\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{n=1}{n\ne m}}^{\tilde{M}-1}{{\rm{\Gamma }}}_{{nm}}\right\}.\end{eqnarray*}$$

Each term in the first summation is given by Equation (8). Because each of the $\tilde{M}$ pairs is correlated with $2(M-2)$ other pairs, the number of nonzero terms in the second sum is $M(M-1)(M-2)$, with each contribution given by Equation (9). The variance then reduces to

$$\begin{eqnarray*}&&{\mathbb{V}}[{\langle Y\rangle }_{\tilde{M}}]=\displaystyle \frac{1+2(M-1){A}^{2}{\cos }^{2}\phi }{\tilde{M}}.\end{eqnarray*}$$

The product ${I}_{a}{I}_{b}$ has similar averaged properties but with the cosine terms replaced by sines. Since all sine and cosine terms are squared, it is easy to see how they drop out when the two averaged products are added to a single average co-spectrum.

Finally, considering an ensemble of L separate observations, the co-spectrum estimator gives a mean and variance of

$$\begin{eqnarray}&&{\mu }_{C}={A}^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\sigma }_{C}^{2}=\displaystyle \frac{4+4(M-1){A}^{2}}{{LM}(M-1)}.\end{eqnarray} \tag{ 11 }$$

Given a sample of K observations with a two-detector system, I can distinguish three averaging scenarios. First, averaging in the time domain gives a standard deviation from Equation (6) using M = 2 and L = K. Second, averaging by correlating in the Fourier domain gives a standard deviation from Equation (11) using M = 2 and L = K. Third, averaging powers gives a standard deviation as in Equation (6) with M = 1 and $L=2K$. For pure noise, the ratios of these uncertainties relate as 1:$\sqrt{2}$:$\sqrt{2}$, indicating that time-domain averaging gives the best detection sensitivity, whereas correlating the concurrent time series gives the same sensitivity as combining the power spectra.