SPIN-DOWN OF THE LONG-PERIOD ACCRETING PULSAR 4U 2206+54

First observed by Uhuru and Ariel V, the low-luminosity (∼10³⁵ erg cm⁻² s⁻¹) hard X-ray source 4U 2206+54 is one of a handful of well-studied galactic high-mass X-ray binaries for which the nature of the accreting compact star is uncertain.

Optical and UV spectroscopy have showed the optical counterpart, BD +53° 2790, to be an unusual O9 active star with strong wind resonance lines in the ultraviolet (Negueruela & Reig 2001). Several authors have concluded from the high column depth and the flaring observed in the X-rays that the compact object is accreting from a strong stellar wind. Ribó et al. (2006) determine from the IUE spectra a relatively low wind velocity of ∼350 km s⁻¹.

Analysis of the first five years of the RXTE ASM light curve of 4U 2206+54 revealed a periodic modulation at a period of 9.57 days, which was interpreted as the orbital period of the binary system (Corbet & Peele 2001). Radial velocity measurements were made for BD +53° 2790; however, it was not possible to independently determine the orbital period from them (Blay 2005). Recent analyses of the Swift/BAT and RXTE/ASM light curves have thrown this picture into confusion (Corbet et al. 2007). The recent data show a modulation with a period of 19.25 days—twice what was previously thought to be the period. If 19.25 days is the orbital period, then the profile of the orbital modulation has evolved from having two peaks per orbit to one peak per orbit.

Using observations extended over 7 days from the Proportional Counter Array (PCA) on the Rossi X-ray Timing Explorer (RXTE), Reig et al. (2009) discovered pulsations at a period of 5559 ± 3 s. Detecting pulsations in this period range from a low Earth orbit is usually difficult because the pulse period is near the orbital period, resulting in poor sampling of pulses. Their data however include several intervals which were uninterrupted by Earth occultation. They also confirmed this discovery using observations from IBIS on the International Gamma-ray Astrophysics Laboratory (INTEGRAL) which has a high Earth orbit with a three-day period.

In Section 2, we present results from new observations of 4U 2206+54 that we have made with Suzaku. These observations were contemporaneous with those of RXTE (Reig et al. 2009). We find strong pulsations at a pulse period of 5554 ± 9 s, in agreement with the RXTE results. In Section 3, we present a reanalysis of the BeppoSAX observations of 4U 2206+54 which show strong pulsations at a period of 5420 ±  28 s; in Section 4, we discuss simulations of Suzaku and BeppoSAX data, in Section 5 the EXOSAT observations of 4U 2206+54, and in Section 6 we briefly analyze the evolutionary status of the system and the nature of its primary component. We conclude that the primary component of the system is a neutron star, which has a huge (∼(3–5) × 10¹⁵ G) surface field and is accreting material onto its surface from a disk. In this light, 4U 2206+54 appears to be the first accretion-powered magnetar identified so far.

We obtained an observation of 4U 2206+54 with Suzaku with the primary objective of searching for long-period pulsations with the X-ray Imaging Spectrometer (XIS; Komaya et al. 2007). To minimize the number and size of data gaps, the observation was time constrained to a phase of the spacecraft precession cycle where the source was always above the Earth horizon. The observation occurred 2007 May 16–17 (MJD 54236.176−54237.813), lasting for 141 ks, with a total exposure for the XIS of 114.8 ks after filtering. A 20.9 ks duration portion of this observation (MJD 54236.861−54236.110) overlaps with the RXTE observations presented by Reig et al. (2009).

The XIS observations were taken in the 1/4 window mode to avoid potential problems with pile-up. Our analysis began with the version 2.0.6.13 processing data (which have aspect correction for thermal wobbling). The charge transfer inefficiency (CTI) correction was remade, and the event data cleaned using the standard filtering with the exception that the lower bound on the dayside Earth limb elevation was reduced from 20° to 175 (which did not appear to increase the broadband background). Events were selected from the active telescopes (XIS0, XIS1, and XIS3) using 208'' × 261'' rectangular extraction regions centered on the source, with the long axis aligned with the 1/4 window. As a limit on the background, events were also extracted from each telescope in regions of the same area as distant as possible from the source. The average rate summed over the telescopes was 13.9 counts s⁻¹ for the source regions and 0.125 counts s⁻¹ for the background regions.

The XIS light curve in the 0.5–10.0 keV range with ∼200 s binning is shown in Figure 1. In order to avoid the loss of data, unequal bin sizes were used. Each good-time interval was divided into equal width bins, with a number of bins chosen to result in the bin width closest to 200 s. There appears to be pulsation with a period near 0.06 days, which is most evident in the spacing of the minima. Also evident is strong variability at much longer timescales.

Zoom In Zoom Out Reset image size

2.1. Pulse Search Statistic

To search for pulsation and determine its period, we have developed a test statistic, which like the Lomb test statistic is based on hypothesis testing, but which can account for the low-frequency variablity and other characteristics of the observations.

The Lomb method (Lomb 1975), a standard technique for detection of pulsations in unevenly sampled data, assumes a null hypothesis of Gaussian white noise, and an alternate hypothesis of Gaussian white noise plus a sinusoid. For 4U 2206+54, the assumption of Gaussian white noise is inappropriate. Reig et al. (2009) show that the power spectra of the flux in the 2–20 keV band are power laws above 1 mHz, with an index between 1.5 and 1.9. The light curve in Figure 1 shows both short-term flaring and strong long-period (∼day) variability, suggesting these red noise spectra extend to lower frequencies. In addition to the power spectrum not being white, the distribution of the rates is far from being Gaussian. In the top panel of Figure 2, we show the distribution of the rates measured with Suzaku XIS, which is more exponential than Gaussian. In the bottom panel, we show the distribution of the natural log of the rate, which has a more symmetric distribution which is closer to a Gaussian.

Zoom In Zoom Out Reset image size

For a more realistic null hypothesis we will assume that the natural log of the flux is Gaussian noise, with a power spectrum of the form

$$\begin{equation} \mathcal {S}(f) = s_0 f_0^\Gamma (f^2+\alpha ^2)^{-\Gamma /2}, \end{equation} \tag{ 1 }$$

where s₀, α, and Γ are adjustable constants and f₀ is a fixed reference frequency. The covariance of the natural logs of rates k and l is then

$$\begin{equation} V_{kl} = \sigma _k^2\delta _{kl}+\int _0^\infty \mathcal {S}(f) {\rm sinc}^2(\pi f \tau) \cos (2\pi f [t_k-t_l]) df, \end{equation} \tag{ 2 }$$

where σ_k is the error on y_k = ln(r_k × 1s), τ is the average rate bin width, and the sinc² term accounts for the effect of binning. The likelihood of the null hypothesis is then given by

$$\begin{equation} \mathcal {L}_0 = {\rm Det}(2\pi V)^{-\frac{1}{2}} \exp \left(-\frac{1}{2}\sum _{kl} [y_k-\mu ] V^{-1}_{kl} [y_l-\mu ]\right), \end{equation} \tag{ 3 }$$

where μ is the mean natural log rate. In practice, it is simpler for us to use the Cash statistic (Cash 1979) $\mathcal {C} =$ $-2{\rm ln}\mathcal {L}$, giving

$$\begin{equation} \mathcal {C}_0 = \sum _{kl} (y_k-\mu) V^{-1}_{kl} (y_l-\mu)-{\rm ln}({\rm Det}(V))+\mbox{const}, \end{equation} \tag{ 4 }$$

which is to be minimized with respect to s₀, a₀, Γ, and μ.

For the alternative hypothesis, the natural log of the flux is Gaussian noise with power spectra of the same form plus a sinusoidal profile. The resulting Cash statistic is

$$\begin{equation} \mathcal {C}_1 = \sum _{kl} (y_k-p_k) V^{-1}_{kl} (y_l-p_l)-{\rm ln}({\rm Det}(V))+\mbox{const}, \end{equation} \tag{ 5 }$$

where the pulse profile p is given by

$$\begin{equation} p_k = \mu + a\cos {2\pi f t_k}+ b\sin {2\pi f t_k}. \end{equation} \tag{ 6 }$$

$\mathcal {C}_1$ is to be minimized with respect to s₀, a₀, Γ, μ, a, and b for a given trial frequency f. The test statistic analogous to the Δχ² is

$$\begin{equation} \Delta \mathcal {C} = \min (\mathcal {C}_0)-\min (\mathcal {C}_1). \end{equation} \tag{ 7 }$$

In the Appendices, we explain how the calculation of this statistic is implemented numerically. In Section 4, we describe how results from this statistic compare with the Lomb method.

2.2. Suzaku Timing Results

In Figure 3, we show $\Delta \mathcal {C}$ for the Suzaku XIS light curve. The false detection probability of the highest peak is 5 × 10⁻¹² including the number of search trials within the 0–1.0 mHz interval.

Zoom In Zoom Out Reset image size

We have investigated the presence of harmonics of the pulse frequency by including additional sinusoids in the profile (Equation (6)). Adding first harmonic terms decreases the Cash $\mathcal {C}$ statistic by 9.9, which has a chance probability of 0.7%. Adding second harmonic terms decreases the Cash statistic by an additional 11.4, which has a chance probability of 0.3%. Further harmonics continue to decrease the Cash statistic, but these decreases are not individually as significant. A period of 5554 ±  9 s and power spectral parameters determined with these two harmonics included are listed in Table 1.

Table 1. Period and Power Spectra Parameter Estimates

Observation	Mid-time	Period	S0a	α	Γ	Harmonicsb
	(MJD)	(s)	(Hz−1)	(10−5 Hz)
Suzaku XIS	54237.0	5554 ± 9	31.6 ± 1.9	2.2 ± 1.5c	1.10 ± 0.08	3
BeppoSAX MECS	51141.1	5420 ± 28	21.4 ± 3.6	<12d	1.14 ± 0.18	2

Notes. ^aPower spectra normalization at frequency f₀ = 1 mHz; ^bNumber of harmonics in fit including fundamental; ^cThe 90% confidence range is 1.6 × 10⁻⁶–5.1 × 10⁻⁵ Hz; ^d90% confidence upper limit.

Download table as:  ASCIITypeset image

In Figure 4, we show the epoch folded Suzaku XIS rates using the period determined from the maximum likelihood fit with two harmonics included. The errors used are the square root of the sample variance in each phase bin divided by the number of rates in the bin. Using the pulse period, power spectrum, and profile parameters estimated from this same fit, we made 10,000 simulated light curves with the same binning as in the Suzaku data. The simulation method is discussed in the Appendices. The smooth curve in Figure 4 is the average epoch folded profile of these simulated light curves.

Zoom In Zoom Out Reset image size

In Figure 5, we show pulse profiles obtained from the Suzaku XIS and HXD/PIN data. Background subtracted light curves with 24 s resolution were made for each energy band and epoch folded with a period of 5554 s. For the PIN profiles, light curves were made from the cleaned event files, with live-time correction based on 144 s averages of the pseudoevents. The subtracted background was from the version 2 non-X-ray background file plus the cosmic X-ray background calculated with the flat response per the recipe on the Suzaku Data Analysis web page.⁶ For both XIS and PIN profiles, the errors were calculated from the sample variance in each bin. Short flares in the data have resulted in the errors being correlated between neighboring bins.

Zoom In Zoom Out Reset image size

The modulation fraction of the profiles ([max-min]/[max+ min]) is 49% in the 0.5–2.5 keV energy band, increasing to 55% in the 20–30 keV band. A notch at the minimum is apparent in all profiles up to 30 keV. The profile is symmetric at low energies, with the leading edge becoming stronger than the tailing edge at higher energies.

We show in Figure 6 how the estimate of the pulse profile and power spectrum of the natural log of the rates is related to the empirically determined power spectrum of the rates at higher frequencies. For the higher-frequency power spectrum, we used a 37 ks segment of the Suzaku XIS observations which contained no gaps. We made a light curve of this segment with 2 s resolution and constructed from it power spectrum estimates based on the Fourier transform of the counts. We then produced 10,000 simulated light curves using the parameters estimated in the maximum likelihood fit of the Suzaku XIS data. Unlike the previous simulations, these light curves had uniform sampling covering the duration of the observation. The average power spectrum from these simulations is shown by the solid curve in Figure 6.

Zoom In Zoom Out Reset image size

We have reanalyzed the BeppoSAX observation presented by Masetti et al. (2004). The source was observed with BeppoSAX on 1998 November 11. We extracted events in the 1.65–10.0 keV range from the version 2 merged Medium Energy Concentrator Spectrometer (MECS) 2 and 3 files obtained from the ASDC Multi-Mission Interactive Archive, using a 4' radius extraction circle centered on the source. As a limit on the background, events were also extracted from an annulus with the same area centered on a radius of 6'. The mean rate from the source region was 0.41 counts s⁻¹ while that of the background region was 0.014 counts s⁻¹, with much of this likely associated with the wings of the point-spread function.

In Figure 7, we show the MECS light curve with ∼200 s binning of the 1.65–10 keV events. Again there appears to be pulsation with a period near 0.06 days, which is most evident in the spacing of the minima.

Zoom In Zoom Out Reset image size

We show a search for long-period pulsations in Figure 8 using the $\Delta \mathcal {C}$ statistic with the same modeling used for the Suzaku observation. The highest peak has a probability of occurring by chance of 5 × 10⁻⁸ including the number of search trials in the 0–1.0 mHz interval. From the location of the highest peak, we estimate a pulse period of 5393 ± 28 s. The peak near 0.36 mHz suggests a first harmonic. We investigated the presence of harmonics of the pulse frequency by including additional sinusoids in the profile (Equation (6)). Adding first harmonic terms to the pulse profile decreases the Cash statistic by 8.3 (1.6% chance probability). Improvements from adding additional harmonics were not as significant. A period estimate of 5420 ± 25 s is obtained from the fit including the first harmonic. This and the estimated power spectrum parameters are shown in Table 1.

Zoom In Zoom Out Reset image size

We find small but sometimes significant differences between the periods estimated with our maximum likelihood analysis and those determined with the Lomb method. We also find differences between the period determined with the maximum likelihood method when harmonics are included with the fundamental in the fit and the period determined when only the fundamental is included. Here we present simulations which examine the extent to which these differences are expected. Our simulation method is presented in the Appendices.

In Figure 9, we present the results of a search for pulsation in the XIS light curve using the Lomb method (Lomb 1975). For the Lomb power we are using twice the normalization given in Press et al. (1992, p. 569). A significant peak is evident near 0.18 mHz. Using the assumptions standard for the Lomb periodigram, the chance probability of this peak to occur within the 0–1.0 mHz range due to noise is 2 × 10⁻³¹. However, these assumptions are not correct as we have discussed, and this significance is overestimated. This should be compared with Figure 3. Note the decrease of the power of the lowest frequency peaks relative to the main peak using the ΔC statistic.

Zoom In Zoom Out Reset image size

The pulse period estimated from the peak of the Lomb power is 5541 ± 11 s which is nearly consistent with our result of 5554 ± 9 s using the maximum likelihood fit with the first and second harmonics included. We simulated 10,000 Suzaku light curves using the same times and binning as the light curve in Figure 1, and the pulse period, power spectra, and Fourier coefficients obtained from our maximum likelihood fit for this light curve. For each, we used the Lomb method to estimate the pulse period. The period distribution was consistent with a Gaussian with mean displaced by −1.8 ± 2 s from the simulated period of 5554 s, and a standard deviation of 17.6 ± 1 s. So for the Suzaku observations we conclude that the Lomb method is essentially unbiased, but with underestimated errors. The period measured with the maximum likelihood fit including only the fundamental was consistent within error to our final result from the fit including the first and second harmonics.

To perform the same investigation for BeppoSAX data, we made simulated light curves with the time binning of the MECS light curve using the estimated parameters from the maximum likelihood fit including the first harmonic. Using 10,000 simulated light curves we found that the distribution of the pulse period determined with the Lomb method was Gaussian with a mean displaced −23.8 ± 0.3 s from the simulated period (5420 s) and a standard deviation of 32 s. The period obtained using the Lomb method and using the actual data is displaced by −44 s from the period from our final maximum likelihood fit, consistent with the distribution observed in the simulations.

The period determined with the maximum likelihood fit of the BeppoSAX data using only the fundamental was displaced by −27 s from the period of 5420 ± 25 s determined from our final fit which included the first harmonic. Using 200 of the simulated light curves, we found the distribution of the periods determined with our maximum likelihood method with only the pulse fundamental had a mean displaced from the simulated period (5420 s) by −17 ± 2 s and a standard deviation of 28 ±  2 s.

Using 200 simulated BeppoSAX light curves, we found no bias in the periods estimated with our maximum likelihood method with the fundamental and first harmonic included. However, the standard deviation of the estimated periods was 28.4 ± 1.4 s which is larger than the error (25 s) we obtained from the curvature of the Cash statistic versus frequency curve for the actual data. This is reflected in the error given in Table 1.

Observations of 4U 2206+54 were made with the medium energy (ME) proportional counters on the European Space Agency's X-ray Observatory EXOSAT. EXOSAT/ME observed 4U 2206+54 on three occasions: 1993 August 8, 1984 December 7, and 1985 June 27. These observations were originally presented by Saraswat & Apparao (1992).

Reig et al. (2009) analyzed the ME standard product light curves from these observations, which are available at HEASARC, and reported pulsation with a 5525 ± 30 period. In Figure 10, we show the light curves for the first and third observation with 200 s binning. We concluded that the background was not correctly subtracted for the standard products light curve of the second observation, and therefore have not shown it.

Zoom In Zoom Out Reset image size

Given the short durations of these light curves (9630 s, 9400 s, and 7404 s, respectively) relative to the pulse period, we do not believe a period measurement based on them is reliable. To demonstrate this, we have made simulations based on two long gapless intervals of the Suzaku XIS light curve. For the first simulation we selected a number of overlapping 9600 s data segments within each gapless interval and from each segment made a pulse period estimate using the Suzaku XIS rates in that segment. Since Reig et al. (2009) used the CLEAN method for period determination we did also. CLEAN attempts to iteratively reconstruct the Fourier amplitude spectrum of the data (Roberts et al. 1987). We used CLEAN within the program PERIOD version 5.0–2, distributed with the Starlink Software Collection. In our analysis of each segment, we used five iterations with a loop gain of 0.2. In the top panel of Figure 11, we show the period from each segment versus the midpoint time of the segment. The results of the second simulation, which used 7400 s segments, are shown in the bottom panel. It is clear that the period estimates from these short data segments vary widely from the 5554 s period determined from the whole observations, and that the formal errors, which do not account for the underlying red noise, grossly underestimate the error distribution.

Zoom In Zoom Out Reset image size

We interpret the observed X-ray pulsations of 4U 2206+54 as due to the rotation of an accreting magnetized neutron star. An alternative interpretation, which was also discussed for the 10,000 s pulsations of the low-luminosity wind-fed system 2S 0114+650 (Corbet et al. 1999), is that the observed period originates from a pulsation of the companion BD +53° 2790. However, this explanation seems less plausible because such persistent modulations should be detectable in the optical photometry, but have not been reported. A mechanism for transferring the modulations to the stellar wind that results in a large modulation of the X-ray flux has not been advanced. In addition, the evolution of the pulse profile with energy that we see in the Suzaku data would not be expected.

Is it reasonable for us to assume that such a long-period accreting pulsar can be formed? The longest period which a neutron star can reach in its evolution is the period at which it emerges from the subsonic propeller state (Ikhsanov 2007):

$$\begin{equation} p_{\rm br} = 15000\ \mu _{32}^{16/21} m^{-4/21} \left(\frac{\dot{M}_{\rm c}}{10^{15}\;{\rm g\;s^{-1}}}\right)^{-5/7}\;{\rm s}, \end{equation} \tag{ 8 }$$

where μ₃₂ and m are the dipole magnetic moment and mass of the neutron star in units of 10³² G cm³ and 1.5 M_☉, respectively. $\dot{M}_{\rm c} = \pi R_{\rm G}^2 \rho V_{\rm rel}$ is the mass with which a neutron star of mass M_ns interacts in a unit time moving through the wind of an average density ρ with a relative velocity V_rel, and R_G = 2GM_ns/V²_rel is the Bondi radius. For the Suzaku observations, we find a 0.5–70 keV luminosity of 4 × 10³⁵ erg s^-1 implying a mass accretion rate onto the stellar surface of $\dot{M}_{\rm a} = r_{\rm ns} L_{x}/GM_{\rm ns} \simeq 2 \times 10^{15}\;{\rm g\;s^{-1}}$, assuming a distance of 2.6 kpc (Blay et al. 2006). The mass accretion rates we inferred from published observations all fall in the range of 4 × 10¹³–3 × 10¹⁵ g s^-1, with the average near the higher end. Assuming $\dot{M}_{\rm c} \simeq \dot{M}_{\rm a}$ (i.e., direct accretion scenario) and putting this into Equation (8), one finds that the condition p_br ≳ 5500 s is satisfied only if the surface field of the neutron star in the present epoch is B₀ ≳ 10¹⁴ G.

An independent estimate of the field strength can be found considering the neutron star spin evolution. The average rate of the spin frequency change between the BeppoSAX and Suzaku observations was $\dot{\nu } = (-1.7 \pm 0.3) \times 10^{-14}\;{\rm Hz\;s^{-1}}$. This indicates that the average spin-down torque applied to the star during this time was $K_{\rm sd} \gtrsim 2 \pi I |\dot{\nu }|$, where I is the star's moment of inertia. If we now adopt the canonical prescription for the spin-down torque (K_sd = k_tμ²/r³_c; see, e.g., Lynden-Bell & Pringle 1974; Lipunov 1992), one finds μ ≳ μ_m, where

$$\begin{equation} \mu _{\rm m} \simeq 10^{32}\ k_{\rm t}^{-1/2} m^{1/2} I_{45}^{1/2} \dot{\nu }_{-14}^{1/2} \left(\frac{p}{5500\;{\rm s}}\right)\;{\rm G\;cm^3}. \end{equation} \tag{ 9 }$$

Here $\dot{\nu }_{-14} = |\dot{\nu }|/10^{-14}\;{\rm Hz\;s^{-1}}$ and I₄₅ = 10⁴⁵ g cm². r_c = (GMp²/4π²)^1/3 is the corotation radius of the neutron star and k_t is a dimensionless parameter limited to k_t ≲ 1. The value of μ_m represents the lower limit to the dipole magnetic moment of the neutron star since the spin-up torque in the above calculations was assumed to be negligibly small. Therefore, our estimate remains valid independently of whether the star between the BeppoSAX and Suzaku observations was persistently in the accretor state or its state was temporarily changed to the propeller. Thus, the 5500 s pulsations in the X-ray flux of 4U 2206+54 can be explained in terms of the spin period of the degenerate companion provided the neutron star in this system is a magnetar whose surface field at the present epoch exceeds 10¹⁴ G.

The modeling of 4U 2206+54 in terms of an accretion-powered magnetar suggests that the neutron star is relatively young (its age is limited to the characteristic time of the supercritical magnetic field decay), with an extended magnetosphere (r_m ∝ μ^4/7), and a relatively small area of hot polar caps at the base of the accretion column (A_p ∝ μ^−8/7). However, the assumption about spherical geometry of the accretion flow, used for making estimates (8) and (9), encounters major difficulties in explaining the mode by which the accretion flow enters the magnetic field of the neutron star at the magnetospheric boundary. As first shown by Arons & Lea (1976) and Elsner & Lamb (1977), a steady accretion onto the stellar surface in this case could occur only if the X-ray luminosity of the pulsar meets the condition L_x>L_cr, where

$$\begin{equation} L_{\rm cr} \simeq 10^{37}\ \mu _{32}^{1/4}\ m^{1/2}\ r_6^{-1/8}\;{\rm erg\;s^{-1}}. \end{equation} \tag{ 10 }$$

Here r₆ is the radius of the neutron star in units of 10⁶ cm. Otherwise, the interchange instabilities of the magnetospheric boundary are suppressed, and the entry rate of accretion flow into the magnetosphere is too small to explain the X-ray luminosity of the pulsar. The star in this case would appear rather as a burster than a persistent source (Lamb et al. 1977).

It is easy to see, however, that 4U 2206+54 does not meet the above criterion. Nevertheless, it is a persistent accretion-powered X-ray pulsar with a relatively small amplitude of variations in X-ray intensity. This contradiction may indicate that either the spherical accretion is realized in a settling (Narayan & Medvedev 2003; Ikhsanov 2005) rather than the direct accretion mode or the neutron star is accreting material from a disk.

A scenario in which the neutron star accretes material from a hot envelope in a settling mode on the bremsstrahlung cooling timescale implies the following condition to be satisfied (Ikhsanov 2005):

$$\begin{equation} \dot{M}_{\rm c} \gtrsim \left[\frac{t_{\rm br}(r_{\rm m})}{t_{\rm ff}(r_{\rm m})}\right]\ \frac{L_{x} r_{\rm ms}}{GM_{\rm ns}}. \end{equation} \tag{ 11 }$$

Here t_ff(r_m) and t_br(r_m) are the free-fall time and bremsstrahlung cooling time at the magnetospheric boundary. Using parameters of 4U 2206+54, one finds $\dot{M}_{\rm c} \gtrsim 10^{17}\;{\rm g\;s^{-1}}$. Therefore, the persistent behavior of the source can be explained in terms of the settling accretion scenario only if the strength of the wind overflowing the neutron star exceeds the typical mass-transfer rate in a binary system with a O9.5 V by almost four orders of magnitude (Reig et al. 2009). Currently available observations do not favor the assumption about so intensive outflow from the normal companion. However, if it were valid, the condition p_br>5500 s could be satisfied only for the surface field of the neutron star ≳4 × 10¹⁵ G.

It is interesting that a similar value of the magnetic field can be found considering a situation in which the neutron star accretes material from a disk. A spin-down of the neutron star in this case indicates that its spin period is smaller than the equilibrium period, p_eq, which is defined by equating the acceleration, $K_{\rm su} \,{=}\, \dot{M}\sqrt{GM_{\rm ns}r_{\rm m}}$, and deceleration, K_sd = k_tμ²/r³_c, torques. Using parameters of 4U 2206+54, one finds that the condition p ≲ p_eq is satisfied only if the strength of the dipole field at the stellar surface in the present epoch is limited to

$$\begin{equation} B(r_{\rm ns}) \gtrsim 3 \times 10^{15}\ \kappa _{0.5}^{7/24} k_{\rm t}^{-7/12} m^{5/6} \dot{M}_{15.3}^{1/2} \left(\frac{p}{5500\;{\rm s}}\right)^{7/6}\ {\rm G}, \end{equation} \tag{ 12 }$$

where κ_0.5 = κ/0.5 is the parameter accounting for the geometry of the accretion flow normalized following Ghosh & Lamb (1978). The spin-down timescale of a newly formed neutron star to a period of 5500 s under these conditions is close to 5000 yr (see Equation (20) in Ikhsanov 2007). This exceeds the decay time of the magnetic field confined to the crust of the neutron star (Colpi et al. 2000), but is still smaller than the decay time of the field permeating the core. Since the X-ray spectrum of the source resembles spectra of accretion-powered pulsars, one can assume that the contribution of the source associated with the field decay to the system X-ray luminosity is small. The field decay time in this case can be limited to

$$\begin{equation} t_{\rm dec} \gtrsim 2 \times 10^5\ \mu _{33.3}^2\ r_6^{-3}\ L_{35.6}^{-1}\;{\rm yr}, \end{equation} \tag{ 13 }$$

which is consistent within results of Heyl & Kulkarni (1998), who treated the magnetothermal evolution of the star in terms of ambipolar diffusion. Here L_35.6 is the X-ray luminosity of the source expressed in units of 4 × 10³⁵ erg s^-1. The proton cyclotron feature in this case is expected to be observed at the energy of ∼30(B₀/5 × 10¹⁵ G) keV.

It, therefore, appears that our interpretation of the observed 5500 s pulsations in terms of spin period of the neutron star is reasonable and suggests that the neutron star in 4U 2206+54 is a magnetar with a surface field of about (3–5) × 10¹⁵ G. A question about the geometry of the accretion flow remains, however, open so far. The spherical approximation of the flow geometry can be applied only under the condition $\dot{M}_{\rm c} \gg \dot{M}_{\rm a}$. On the other hand, the observed wind velocity of the O-star companion, V_w ≳ 350 km s^-1 (see Ribó et al. 2006), substantially exceeds the upper limit to the relative velocity between the star and the wind at which the disk formation could be expected:

$$\begin{equation} V_{\rm rel} \lesssim V_{\rm cr} \simeq 120\ \xi _{0.2}^{1/4}\ \mu _{33}^{-1/14}\ m^{11/28}\ \dot{M}_{15}^{1/28}\; P_{20}^{-1/4}\;{\rm km\;s^{-1}}. \end{equation} \tag{ 14 }$$

Here P₂₀ is the orbital period of the system in units of 20 days, and ξ_0.2 = ξ/0.2 is the parameter accounting for the inhomogeneities of the accretion flow, which is normalized following Ruffert (1999). Nevertheless, one cannot exclude the presence of a fossil disk surrounding the magnetosphere of the neutron star which could be formed at the fall-back stage of the magnetar evolutionary track (Woosley 1988; Ekşi & Alpar 2003). Further studies of appearance of such a disk in a situation when the neutron star rotates extremely slowly and has a huge magnetosphere could help to choose appropriate flow geometry and reconstruct the accretion picture in this enigmatic source.

Finally, we would like to point out here that at the rate $\dot{\nu } = (-1.7 \pm 0.3) \times 10^{-14}\;{\rm Hz\;s^{-1}}$ the pulsar would reach zero frequency in 300 yr, making it unlikely that this is a long-term trend. Over long timescales, wind-fed pulsars tend to show random walk behavior in their spin frequency. For 4U 2206+54 we can estimate the random walk strength as S = <(Δν)²/Δt> ≃ 8 × 10⁻²⁰ Hz² s^-1. For comparison, Vela X-1 has a random walk strength of 2 × 10⁻²⁰ (Deeter et al. 1987). The random walk behavior must break down, and the torque becomes correlated at short timescales since the spin-up rate $\dot{\nu } \sim (S/\Delta t)^{1/2}$ from wind accretion is limited by that due to disk accretion at the same mass accretion rate. Assuming a mean mass accretion rate of 10¹⁵ g s^-1, and a magnetic field of 10¹⁵ G, we find the spin-up rate must be correlated on timescales shorter than ∼5 days.

Based on the strong emission in helium lines, Negueruela & Reig (2001) and Blay et al. (2006) have proposed that the optical companion BD +53° 2790 is similar to θ¹ Ori C, which has a oblique magnetic dipole with kG surface strength which magnetically channels its stellar wind. If so, this long correlation timescale could be related to long-term correlations in the wind speed or density associated with this magnetic channeling.

M.H.F. acknowledges support from NASA grant NNX08AG12G. N.R.I. acknowledges support from NASA Postdoctoral Program at NASA Marshall Space Flight Center, administered by Oak Ridge Associated Universities through a contract with NASA, and support from Russian Foundation of Basic Research under the grant 07-02-00535a.

Numerically, we calculate V in Equation (2) as

$$\begin{equation} V_{kl} = \sigma _k^2\delta _{kl} +\sum _{j=0}^M w_j \mathcal {S}(f_j) {\rm sinc}^2(\pi f_j \tau) \cos (2\pi f_j [t_k-t_l]) \Delta f. \end{equation} \tag{ A1 }$$

Here the frequency f_j = jΔf and the frequency step used is Δf = (16T)⁻¹, with T being the duration of the observation. The weights w_j give trapezoidal integration with w_j = 1 except for w₀ = w_N = 0.5. M is chosen as the nearest integer to T/τ so that the integral extends to where the power is cut off by the binning effect.

By using

$$\begin{eqnarray} \cos (2\pi f_j [t_k-t_l]) &=& \cos 2\pi f_j t_k \cos 2\pi f_j t_l \nonumber \\ && +\sin 2\pi f_j t_k \sin 2\pi f_j t_l, \end{eqnarray} \tag{ A2 }$$

Equation (A1) can be expressed as

$$\begin{eqnarray} &&V = U^T U, \end{eqnarray} \tag{ A3 }$$

where the matrix U has N + 2M + 1 rows and N columns, with N being the number of measurements, and T signifies the matrix transpose. Explicitly, U is given by

$$\begin{equation} \begin{array}{@{}ll} \qquad U_{i,k} = \sigma _i \delta _{i,k} \; {\rm for }\; 0 \le i < N\\ \quad \,\,\, \,\,\,U_{N,k} ={[w_0 \mathcal {S}(0) \Delta f]^\frac{1}{2}}\\ U_{N+2j-1,k} ={[w_j \mathcal {S}(f_j) {\rm sinc}^2(\pi f_j \tau)\Delta f]}^\frac{1}{2}\cos 2\pi f_j t_k\\ \qquad \qquad \quad {\rm for }\ 1 \le j \le M\\ \quad U_{N+2j,k}={[w_j \mathcal {S}(f_j) {\rm sinc}^2(\pi f_j \tau)\Delta f]}^\frac{1}{2}\sin 2\pi f_j t_k\\ \qquad \qquad \quad {\rm for }\ 1 \le j \le M. \end{array} \end{equation} \tag{ A4 }$$

By applying a series of N orthogonal Householder transformations (reflections), U can be transformed into an N × N upper triangular matrix $\tilde{U}$, while preserving Equation (A3) (Householder 1958). We then have

$$\begin{equation} \mathcal {C}_1 = |\tilde{U}^{-T}(\mathbf {y}-\mathbf {p})|^2 +2\sum _i {\rm ln}(|\tilde{U}_{ii}|), \end{equation}\vspace*{-10pt} \tag{ A5 }$$

with a similar expression for $\mathcal {C}_0$. Here y is the vector with components y_k and p is the vector with components p_k. We will introduce the augmented parameter vector q = [μ, a, b, − 1], and using Equation (6) define the matrix H so that p − y = Hq. Then with $R=\tilde{U}^{-T}H$ we have

$$\begin{equation} \mathcal {C}_1 = |R{\bf q}|^2 +2\sum _i {\rm ln}(|\tilde{U}_{ii}|). \end{equation}\vspace*{-10pt} \tag{ A6 }$$

By applying Householder transformations, R can be transformed into a 4 × 4 upper triangular matrix $\tilde{R}$ while preserving Equation (A6). We can then minimize $\mathcal {C}_1$ with respect to μ, a, and b by setting the first three components of $\tilde{R}{\bf q}$ to zero and solving for μ, a, and b, with the minimum of $|\tilde{R}{\bf q}|^2$ given by R²₃₃.

We minimize $\mathcal {C}_0$ and $\mathcal {C}_1$ with respect to the power spectral parameters s₀, a₀, and Γ by the downhill simplex method.

In our simulations, we first simulate natural log rates y_k with covariance given by Equation (A1) with the measurement errors σ_k set to zero and means p_k given by Equation (6) (perhaps with terms for harmonics). Simulated counting statistical errors are then added to the rates r_k.

As above an upper triangular matrix $\tilde{U}$ is calculated such that $\tilde{U}^T\tilde{U} = V$. Then with x being the vector of N random normal variants with unit variance zero mean, we calculate

$$\begin{equation} {\bPsi } = \tilde{U}^T\mathbf {x}. \end{equation} \tag{ B1 }$$

The covariance matrix of ${\bPsi }$ is

$$\begin{equation} \langle {\bPsi }{\bPsi }^T\rangle = \tilde{U}^T\langle \mathbf {x}\mathbf {x}^T\rangle \tilde{U} = \tilde{U}^T\tilde{U} = V. \end{equation} \tag{ B2 }$$

The simulated rates prior (without counting statistical errors) are given by

$$\begin{equation} {r_k} = \exp ({\Psi _k}+{p_k}). \end{equation} \tag{ B3 }$$

Counting statistical errors are then simulated from these rates and the bin widths.