AN ACCRETION MODEL FOR THE ANOMALOUS X-RAY PULSAR 4U 0142+61

Anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs) are young neutron stars with X-ray luminosities much larger than their spin-down power and long periods in the range 2–12 s. They are widely believed to be magnetars deriving their X-ray emission from the decay of super-strong magnetic dipole fields (B ≳ 10¹⁵ G; e.g., Duncan & Thompson 1992; Thompson & Duncan 1995). At the outset, the magnetar model was developed to explain the giant bursts and the large bursts of SGRs, which exceed the Eddington limit of neutron-star luminosities by a very large factor. If the observed spin-down of some of these sources is interpreted as the consequence of magnetic dipole braking, the resulting polar field strengths are of the order of ≳ 10¹⁵, in qualitative agreement with what has been inferred from the observed luminosities (Kouveliotou et al. 1998). Later on, it was discovered that AXPs show short bursts as well, though less frequent and less energetic than SGRs, leading to the general notion that both types of sources are closely related or represent even a single class. In the magnetar picture, the steady X-ray emission of these sources, which have luminosities of typically a few times 10³⁵ erg s⁻¹, is thought to be caused by a twist of the magnetosphere leading to the amplification of the magnetic field and the acceleration of particles, which produce the X-ray emission. The twist is caused by rotational motions of a crust plate and has a lifetime of ∼1 year (Beloborodov & Thompson 2007). For recent reviews of the magnetar model see Woods & Thompson (2006) and Mereghetti (2008).

An alternative energy source for the persistent and transient X-ray luminosities of AXPs and SGRs is accretion from fallback disks, first proposed by van Paradijs et al. (1995) and followed by Chatterjee et al. (2000) and Alpar (2001). This class of models was developed further in a series of papers (Ekşi & Alpar 2003; Ertan & Alpar 2003; Ertan & Cheng 2004; Ertan et al. 2006, 2007, 2009; Ertan & Çalışkan 2006; Ertan & Erkut 2008). The fallback-disk model gets support from the discovery of IR/optical radiation from two of the AXPs, 4U 0142+61 (Wang et al. 2006) and 1E 2259+586 (Kaplan et al. 2009), which has been successfully interpreted as disk emission. This optical emission is especially important, because the spectral fits at short wavelengths put an upper limit on the inner disk radius and thereby on the magnetic dipole field strength (B < 10¹³ G).

The fallback-disk model explains the spin-down of AXPs and SGRs by the disk–magnetosphere interaction and requires "normal" neutron–star dipole fields (10¹²–10¹³ G). We note that this is in line with the recent discovery of a magnetar (SGR 0418+5729) with a low magnetic dipole field B < 7 × 10¹² G (Rea et al. 2010; see also Alpar et al. 2011). In addition, the model is successful in predicting the period clustering of AXPs/SGRs in the range of 2–12 s. More recently, it was shown that the X-ray luminosity and the rotational properties, including the braking index n = 0.9  ±  0.2, of the so-called high-B radio pulsar PSR J1734−3333 (Espinoza et al. 2011) can be explained by the evolution of a neutron star with a fallback disk and B = 10¹² G (Çalışkan et al., in preparation). On the other hand, this model cannot explain the super-Eddington bursts, which are relatively rare. They are attributed to magnetar-type activities occurring in local multipole fields (star spots).

All references quoted above, regarding the fallback-disk model, are dealing with the physics of the accretion disk and its interaction with the magnetosphere of the neutron star. A first attempt to explain the hard X-ray spectrum of the well-studied AXP 4U 0142+61 by bulk motion and thermal Comptonization has been presented by Trümper et al. (2010). Figure 1 shows the excellent fit of the COMPTB model to the pulse phase-averaged data from Chandra MEG and HEG at low energies and INTEGRAL ISGRI above 15 keV. In the present paper, we study the energy-dependent pulse profiles and the phase-dependent energy spectra, taking again 4U 0142+61 as an example.

Zoom In Zoom Out Reset image size

In Section 2 we summarize the observational properties of 4U 0142+61, in Section 3 we present our model, in Section 4 we interpret the phase-dependent energy spectra and infer the surface magnetic field strength, in Section 5 we discuss the heating of the polar cap, in Section 6 we compare 4U 0142+61 with low-luminosity X-ray pulsars in binary systems, and in Section 7 we give our summary and conclusions.

Since its discovery by Uhuru, 4U 0142+61 has been observed by many groups. The source has a period of 8.69 s and spins down with a rate of 2 × 10⁻¹² s s⁻¹ (Dib et al. 2007). It shows a spectrum in the soft band (0.8–10 keV) that can be fitted by a blackbody (kT = 0.42 keV) and a steep power law (photon-number spectral index Γ_s = 3.3; e.g., Juett et al. 2002). At higher energies (E > 10 keV), a hard spectral tail (Γ_h = 0.93) is observed, which extends beyond 150 keV and contributes about one-third of the total X-ray luminosity (den Hartog et al. 2008). That paper presents the most detailed data on the spectra and light curves of the persistent soft and hard X-ray emission of this source, based on observations with ASCA, XMM-Newton, RXTE, and International Gamma-Ray Astrophysics Laboratory (INTEGRAL). In our present analysis, we use the data provided by den Hartog et al. (2008), since they give very detailed information on the energy-dependent pulse shapes (in particular their Figure 9). According to these data, the luminosity of 4U 0142+61 in the energy band 0.8–160 keV is ∼4.6 × 10³⁵ erg s⁻¹, assuming a distance of 3.6 kpc (Durant & van Kerkwijk 2006). The luminosity in the soft and the hard band is ∼3.2 × 10³⁵ erg s⁻¹ and ∼1.4 × 10³⁵ erg s⁻¹, respectively.

The broadband spectrum of 4U 0142+61 has been measured by Suzaku (Enoto et al. 2011) in the 0.4–70 keV band, confirming the existence of a two-component spectrum. The photon index of the hard component is Γ_h = 0.9. The Suzaku observations have the advantage that data in the hard (>10 keV) and the soft (1–10 keV) bands are taken simultaneously. The luminosities are 2.8 × 10³⁵ erg s⁻¹ in the soft band (1–10 keV) and ∼6.8 × 10³⁴ erg s⁻¹  in the hard band (10–70 keV). We compared the luminosities derived from the Suzaku spectral fits with those from the combined XMM-Newton and INTEGRAL spectral fits in the same bands and found excellent agreement (Table 1).

The same holds for the spectral parameters, except in the case of the hard power-law index derived from Suzaku, which strongly depends on the model for the soft component. In conclusion, this comparison suggests that the quiescent spectrum of the source is quite stable. It also justifies our choice to use the INTEGRAL data, which are of superior statistical quality in the hard spectral range.

Using the observational data of den Hartog et al. (2008), we have derived the observed pulsed profiles by adding the pulsed and non-pulsed components. The result is shown in Figure 2. As noted by den Hartog et al. (2008), morphology changes are ongoing throughout the whole energy range. However, we stress that the principal structure of the pulse profiles is very simple: it mainly consists of two pulses. One is located at phase ∼0.1, extending from the lowest to the highest energies (0.8–160 keV). We call it the main pulse. The secondary pulse is located at phase ∼0.6. Its relative strength decreases strongly from the lowest energy channel (0.8–2 keV) to medium energies (8–16.3 keV) and increases again in the 20–50 keV band.

Zoom In Zoom Out Reset image size

In Figure 3, we show the phase-dependent pulsed spectra derived from the data of den Hartog et al. (2008). The spectrum of the main pulse is taken from the phase interval 0.85–0.35 (shown as a dashed line), while the spectrum of the secondary pulse is taken from the phase interval 0.35–0.85 (shown as a dotted line).

Zoom In Zoom Out Reset image size

The basic picture of our model is shown schematically in Figure 4. It is based on the model that has been developed by many authors, notably Davidson (1973), Basko & Sunyaev (1976, hereafter BS76), and Lyubarskij & Sunyaev (1982, hereafter LS82). The accreting matter falls freely on magnetic field lines until it is stopped in a radiative shock, where it is slowed down from a supersonic to a subsonic flow in the sinking column beneath the shock. Most of the accretion energy released in the radiative shock is transformed into hard X-ray photons by bulk-motion Comptonization (BMC; Blandford & Payne 1981; Payne & Blandford 1981; LS82), which escape sideways, forming a fan beam. We note that in high-luminosity sources, having a large transverse optical depth in the accretion column at the height of the shock, most of the BMC photons will be downscattered before they can escape. As we discuss in Section 3.2 below, due to the low luminosity of AXPs and SGRs, the transverse optical depth is moderate (τ_t ∼ 5). The height of the shock above the stellar surface depends on the mass accretion rate and the cross section of the accretion column. If the fan beam is formed very close to the neutron star surface, then about half of it escapes, while the other half—assisted by gravitational bending—hits the polar cap surrounding the accretion column. The impinging flux heats the polar cap photosphere, and some part of it is scattered off by it. The scattered flux and the photospheric emission form a polar beam, which is orthogonal to the fan beam. If the radiative shock forms at a height above the neutron star surface, then more than half of the fan beam escapes. The larger the height, the larger the fraction of the fan beam that escapes.

Zoom In Zoom Out Reset image size

3.1. Outline of the Model

3.2. Physical Parameters of the Accretion Column

For simplicity, it is usually assumed that the accretion flow is homogeneous and the accretion funnel is cylindrical. However, for detailed comparisons between model and observations of the energy-depended pulse profiles, as we do below, this picture is too simple and qualitatively inadequate.

BS76 and LS82 have pointed out that, for a neutron star with an inclined magnetic field, the cross section of the accretion column at the neutron star surface has the shape of a thin annular arc, located at an off-center position with respect to the magnetic pole. BS76 quote a ratio of width to length of the arc of 0.025. The off-center position of the accretion spot has been confirmed qualitatively by magnetohydrodynamic simulations of Romanova et al. (2004); however, the width of the arc is much larger compared with the idealized model of BS76, namely, ∼0.3–0.5. The simulations of Romanova et al. (2004) have been done for relatively small magnetospheres, like those of accreting millisecond pulsars. Unfortunately, they could not be performed for the large magnetospheres of slowly spinning neutron stars with high magnetic fields due to limitations of computing power (M. M. Romanova 2009, private communication). We assume that the picture of Romanova et al. (2004) is valid qualitatively also in our case. In Figure 4(b), we sketch the footprint of the column at the neutron star surface. This geometry leads to a modulation of the radiating area (and the X-ray flux) with the rotational phase at half the neutron star period. The ratio of minimum to maximum flux is given by d₀/l₀. Due to the opening of the magnetic field lines, the length and the width of the accretion column scale with r as l(r) = l₀ × (r/R)^3/2 and d(r) = d₀ × (r/R)^3/2.

The effective Eddington luminosity in the accretion funnel, near the neutron star surface, is (Becker et al. 2012)

$$\begin{eqnarray} L_{E}^{{\rm eff}} &=& L_{E} \frac{l_0 d_0}{4\pi R^{2}} \frac{\sigma _T}{\sigma _{||}} = \frac{GMm_{p}c}{\sigma _T} \frac{l_0 d_0}{R^2} \frac{\sigma _T}{\sigma _{||}}\nonumber \\ &=& 5.4\times 10^{32} \left(\frac{l_0}{10^{4}{\,\rm cm}}\right) \left(\frac{d_0}{6\times 10^{3}{\,\rm cm}}\right) \left(\frac{M}{1.4\,M_{\odot }}\right)\nonumber \\ && \times \left(\frac{12.5 {\,\rm km}}{R}\right)^{2} \frac{\sigma _T}{\sigma _{||}}\; {\rm erg\,s^{-1}}, \end{eqnarray} \tag{ 1 }$$

where L_E is the Eddington luminosity, M is the mass of the neutron star, σ_T is the Thomson cross section, and σ_|| is the average over energy cross section for photons traveling along the magnetic field.

Since the flow in the accretion funnel is supercritical, a radiative shock is established at a height H above the neutron star surface (BS76, LS82, BW07). The upscattered photons in the shock escape sideways, as a fan beam, and their propagation is controlled by Thomson scattering. The transverse optical depth to electron scattering at or above the shock along the width of the funnel is

$$\begin{eqnarray} \tau _t &=& n_{e}(r) \sigma _T d(r) = {\frac{\dot{M}}{m_p v_{ff}(r) l(r) d(r)} } \sigma _T d(r) \nonumber \\ &=& {\frac{\dot{M}}{m_p v_{ff}(r) l(r) } } \sigma _T ={\frac{\dot{M} \sigma _T}{m_p v_0 l_0 }} {\frac{R}{r}}, \end{eqnarray} \tag{ 2 }$$

where n_e(r) is the electron density, v_ff(r) = v₀(R/r)^1/2 is the free-fall speed, v₀ = (2GM/R)^1/2 ≈ 0.6c(M/M₀)^1/2(12.5km/R)^1/2, and $\dot{M}$ is the mass accretion rate. The transverse optical depth cannot be large, as it is in X-ray pulsars, because the BMC photons would be downscattered before they escaped. On the other hand, it cannot be too small, because BMC would be inefficient. In order to have a transverse optical depth of order a few, for accretion rates of the order of 10¹⁵ g s⁻¹, the length l₀ must be of the order of 10⁴ cm.

The analytic solution for the spectrum produced at the shock, which was found by LS82, provides for a high-energy power-law spectrum with index Γ_h ∼ 0.9. Using a Monte Carlo code (N. D. Kylafis & J. E. Trümper, in preparation), we have obtained a similar spectrum, if the transverse Thomson optical depth at the shock is τ_t(R + H) ∼ 5. Using Equation (2), this implies that for H < R the height of the shock is given approximately by

$$\begin{equation} H \sim \left(\frac{\dot{M} \sigma _T}{5 m_p v_0 l_0} -1\right)R \sim 0.2\,R, \end{equation} \tag{ 3 }$$

where we have used $\dot{M} = 3\times 10^{15}$ g s⁻¹ and l₀ = 10⁴ cm. It is clear from the above that despite the low luminosity of AXPs and SGRs, the radiative shock may be located at some height, depending on the accretion rate, the footprint of the accretion column, and the gravitational field.

For completeness, we mention that in the polar direction (i.e., parallel to the magnetic field), the Thomson optical depth of the accretion column is

$$\begin{equation} \tau _p = \frac{2}{3} \tau _t(R) \frac{R}{d_0} \gg \tau _{t}(R), \end{equation} \tag{ 4 }$$

3.3. Modeling the Pulse Shapes of 4U 0142+61

As already noted above, after almost 40 years of theoretical efforts, there is no fully self-consistent description yet of the radiation from an accretion column. This is true in particular for the angle and energy dependence of the emission. The most advanced approach is the one by BW07, which describes the spectral formation in terms of TC and BMC but uses angle-averaged scattering cross sections, which are not suited to derive angular distributions. We therefore use a pragmatic approach and introduce beaming functions to describe the angular dependence of the different beams. This method has been used widely for the analysis of X-ray pulsar profiles (e.g., Leahy 1990, 2003; Kraus 2001).

To model the observed pulse profiles, we assume that the neutron star has a centered magnetic dipole. Thus, the fan beam, which escapes sideways from the column, and the polar beam are orthogonal. The beaming functions are given by fsin ^mθ and pcos ⁿθ for the fan beam and the polar beam, respectively. Here θ is the polar angle, i.e., the angle between the direction of photon propagation and the magnetic dipole axis. The amplitudes f and p and the exponents m and n are free parameters, which may be different for different energy intervals. Note that m = 1 or n = 1 corresponds to the case of isotropic emission, e.g., blackbody emission. The ansatz for the beaming function is a generalized version of the angular distributions of photons presented by Nagel (1981) and Meszaros & Nagel (1985a, 1985b) for columns and slabs, which can be described by sin ^mθ and cos ⁿθ, respectively, with m and n ⩽ 1.5, except for an energy interval around the cyclotron resonance. In that work, the full angle- and energy-dependent magnetic scattering cross sections, including effects of vacuum polarization, were used, assuming a homogeneous temperature and density distribution in the radiating slab or column. However, since in the outer layers of a slab or a column the plasma temperature will drop toward the surface, photons escaping perpendicular to the surface will escape from deeper layers than photons escaping in other directions. This leads to a limb-darkening effect, which makes the beaming patterns narrower compared with those of Nagel (1981) and Meszaros & Nagel (1985a, 1985b). Thus, the expected exponents of the beaming will be significantly larger than 1.5.

A model with orthogonal fan and polar beams will lead to symmetric pulse profiles with double pulses separated by 0.5 in phase, in disagreement with the observed pulse profiles, which show asymmetries and slightly different pulse separations. To account for these effects, we assume, as discussed in Section 3.2, that the accretion column has a non-cylindrical cross section with major axis l₀ and minor axis d₀ (cf. Figure 4(b)).

To account for gravitational bending, we use the expression derived by Beloborodov (2002). For a photon emitted near the neutron star at an angle θ with respect to the radial direction, the angle θ + δ, at which it arrives at infinity, is given by

$$\begin{equation} 1 - \cos (\theta + \delta) = \frac{1 - \cos \theta }{1-\frac{R_s}{R}}, \end{equation} \tag{ 5 }$$

where R is the radius of the neutron star and R_s = 2GM/c² is its Schwarzschild radius. This equation provides a very good approximation in a Schwarzschild metric for R/R_s > 2. We assume M = 1.4 M_☉ and R = 12.5 km, viz., R/R_s = 3. An important parameter is the height H of the radiative shock above the neutron star surface. In the following, we discuss three cases: H = 0, H = 2 km, and H = 0.5R.

• 1.  
  H = 0: this choice is motivated by the general view that at low accretion rate the matter is stopped by the radiative shock at a position close to the stellar surface. The fitting parameters are the inclination i of the system, the angle α between the rotation axis and the magnetic dipole axis, the axis ratio d₀/l₀, the azimuthal angle ϕ of the long axis, phase zero of the pulse profile, the exponents m and n of the beaming functions for the fan and the polar beam, and the normalizations f and p of the fan beam component and the polar beam component, respectively. We remark that the first five parameters (i, α, d₀/l₀, ϕ, and phase zero) are fitted jointly for all energy bands, while the others are fitted for each energy band separately.The fit is done by a hierarchical grid search. In an outer loop, the global parameters are scanned in a five-dimensional grid. In an inner loop, the two exponents and the flux ratio between the fan and polar beam are scanned in a three-dimensional grid, for each combination of the global parameters of the outer grid. For each candidate pulse profile, the overall normalization of the flux is directly determined by the ratio of the mean flux of that profile and the measured profile in the corresponding energy range.This procedure is applied several times, starting with a coarse grid covering a wide range of parameters. This grid is subsequently refined by a denser grid around the solution found in the preceding step. We use the χ² statistics as the criterion for the agreement between the observed and calculated pulse profiles. In order to avoid that the solution is dominated by the low-energy profiles, which have the best statistical quality, we add systematic errors to the statistical errors, individually for each energy band, so that formally acceptable fits are obtained. This method allows us also to quantify to which accuracy our model is capable of reproducing the observed pulse profiles. We find systematic errors of 2.1%, 1.7%, 4.4%, 3.5%, 3.5% for the pulse profiles of the corresponding ranges 0.8–2.0 keV, 2–4 keV, 4–8 keV, 8–16.3 keV, and 20–50 keV (for the 50–160 keV band, no additional error is required). We then estimate the uncertainty in the model parameters by subsequently varying the individual parameters in both directions (while keeping the other parameters fixed at the best-fit values), until the χ² increases by 1, which corresponds to the 68% confidence level for one interesting parameter. Figure 5(a) shows the model fit to the observed light curves. The inclination is i = 421 ± 04, the angle between the dipole with respect to the spin axis is α = 19° ± 09, and the axis ratio of the accretion column is d₀/l₀ = 0.69 ± 0.01, where the errors are 1σ.Our model gives us information about the full beaming geometry and thus allows us to determine the "intrinsic luminosities" of the different beams. This is a major advantage compared to analysis methods that are purely based on the observed flux, i.e., the flux that we sample along the line of sight during the rotation of the neutron star. Depending on the geometrical circumstances, the observed flux may not be representative for the X-ray emission of the neutron star. In the following, we utilize the information about the full beaming geometry to compute which fraction of the intrinsic fan beam luminosity irradiates the neutron star and compare this fraction with the intrinsic luminosity of the polar beam. For clarity, we focus in the following on the properties of one pole (due to the symmetry of our model, the properties of the other pole are identical).According to the H = 0 model, the intrinsic 0.8–160 keV luminosity of the fan beam is 4.4 × 10³⁵ erg s⁻¹, which includes the fraction that falls onto the neutron star and the one that escapes. The intrinsic luminosity of the polar beam is 0.7 × 10³⁵ erg s⁻¹, viz., 15% of the energy of the intrinsic fan beam. The expected fraction of the total luminosity of the fan beam hitting the neutron star is 50%. Thus, the neutron star reemits only a fraction 15/50 = 30% of the incident flux, while the rest disappears, e.g., in the neutron star crust. As this seems quite unlikely, we investigated an alternative scenario, where the radiative shock is not located at the stellar surface, but at some height H above it. In this case, the neutron star intercepts less than half of the incident fan beam flux.
• 2.  
  H = 2 km: in order to compute which fraction of the total fan beam luminosity reaches the neutron star for the case H > 0, it is necessary to know the trajectories of the emitted photons in the presence of gravitational bending. Such trajectories are found in the work of Leahy (2003) for a neutron star radius R = 12.5 km. He shows that X-rays emitted from the accretion column at a height H reach the neutron star if they are emitted at zenith angles θ > θ_i(H), while the rest escapes. Implementing this condition into the model and assuming H = 2 km, the solution is shown in Figure 5(b). For this model, we find systematic errors of 2.0%, 3.1%, 7.1%, 6.1%, and 7.5% for the pulse profiles of the corresponding ranges 0.8–2.0 keV, 2–4 keV, 4–8 keV, 8–16.3 keV, and 20–50 keV (for the 50–160 keV band, no additional error is required). The inclination is i = 596 ± 24, the angle between the dipole and the spin axis is α = 312 ± 17, and the axis ratio of the accretion column is d₀/l₀ = 4/5. The intrinsic 0.8–160 keV luminosity of the fan beam is 1.8 × 10³⁵ erg s⁻¹, and that of the polar beam is 0.3 × 10³⁵ erg s⁻¹. Thus, the ratio of the luminosities is again 15%. However, the fraction of the total luminosity of the fan beam that reaches the neutron star is now smaller. Its exact value depends on the exponent m. For a typical m = 4, it is about 20%. In this case, an X-ray albedo of 15/20 = 75% would be required for the neutron star, which looks more plausible than the 30% found for H = 0.
• 3.  
  H = 0.5R: we also investigated the case H = 6.25 km and obtained fit results (i = 550, α = 305, d₀/l₀ = 4/5) that are very similar to those for H = 2 km. This is a consequence of the fact that the effect of gravitational bending becomes less important with increasing height. Compared to H = 2 km, the intrinsic luminosity of the fan beam increases slightly (by less than 3%), while the polar beam luminosity stays the same, so that the ratio is again ∼15%. At H = 0.5R, however, the neutron star intercepts only ∼5% of the incident fan beam flux, so that the fan beam could power only 5/15 = 1/3 of the polar beam luminosity.The fact that the fit results are rather insensitive to H, for H between 2 km and 6.25 km, allows us to compute the maximum value of H in a straightforward way, by determining the height at which the neutron star would intercept 15% of the fan beam luminosity. This value depends on the exponent m and decreases with increasing m. For representative values m = 2 and 3, it is ∼3.5 km and ∼2.8 km, respectively. We note that the upper limit of ∼3.5 km or 0.3 R, which we find for the height of the radiative shock, is consistent with the height estimate (∼ 2.5 km) presented in Section 3.2.

Zoom In Zoom Out Reset image size

As already discussed in Sections 3.2 and 3.3 due to the beaming geometry and enhancement by gravitational bending, a substantial part of the fan beam hits the polar cap photosphere and will heat it up. A fraction of the infalling photons will be scattered, leading to a reflected beam emitted in the polar direction. The reflection coefficient depends on photon energy and reaches a maximum at the local cyclotron frequency. Thus, one expects that the intensity ratio of the polar to the fan beam increases toward the cyclotron energy. Such an increase is indeed observed (Figure 3). At ∼60 keV, the polar beam reaches its maximum, followed by a decline or cutoff beyond ∼70 keV. Thus, the corresponding cyclotron energy is ∼60 keV, implying a polar magnetic field of B ∼ 5 × 10¹²(1 + z) G, where z is the gravitational redshift. We note that the rather large apparent width of this feature may be due to the variation of the magnetic dipole field over the polar cap photosphere, which, according to observations, has a radius of ∼8 km (White et al. 1996; Israel et al. 1999; Juett et al. 2002). Clearly, a confirmation of this spectral feature and a more accurate determination of its position and width in energy will allow us to further test the accuracy of these estimates.

It is well known that the polar caps of AXPs/SGRs show significantly higher photospheric temperatures compared with other relatively young neutron stars, like radio pulsars (e.g., Aguilera et al. 2008). It is evident that our model provides a natural explanation for this observational fact, since it predicts a strong illumination of the photosphere surrounding the accretion column by the fan beam. The efficiency of this process is increased by gravitational bending.

Because of the relatively high threshold (>0.8 keV), the XMM-Newton data of den Hartog et al. (2008) taken in 2003 and 2004 do not contain information on the blackbody component. The Chandra MEG and HEG observations of 4U 0142+61 on 2001 May 23 (Juett et al. 2002) have measured the spectrum at lower energies (0.5–10 keV) that can be fitted by a blackbody (kT = 0.418 ± 0.013 keV, L_bb = 5.7 × 10³⁴(d/3.6 kpc)² erg s⁻¹) and a power-law component with a photon index Γ_s = 3.3. The observed energy flux in the 0.8–2 keV interval is 47.5 eV cm⁻² s⁻¹. In our data, based on the XMM-Newton observations, the corresponding flux is 87.6 eV cm⁻² s⁻¹. The difference may be due to source variability and/or systematic errors. Assuming that the ratio of blackbody and power-law energy flux (which is 1) is the same in both data sets, we find a blackbody luminosity L_bb = 10³⁵ erg s⁻¹ expected for the XMM-Newton spectrum.

The blackbody luminosity of 4U 0142+61 can also be derived from data presented by Enoto et al. (2011). They fit two blackbodies to the soft component with kT₁ = 0.337 keV, R₁ = 13.2 km and kT₂ = 0.633 keV, R₂ = 2.3 km. The corresponding blackbody luminosities are L₁ = 7.25 × 10³⁴ erg s⁻¹ and L₂ = 2.69 × 10³⁴ erg s⁻¹. Thus, the total blackbody luminosity is L_bb = 10³⁵ erg s⁻¹, in good agreement with the value derived from the Chandra data for the time of the XMM-Newton observations. This result, together with the agreement between the XMM-Newton data and the Suzaku data presented in Table 1, shows that the observational data are consistent.

Traditionally, the soft component of AXPs and SGRs has been fitted by two different models: (1) blackbody plus power law (BB+PL) or (2) double blackbody (dBB). Halpern & Gotthelf (2005) discussed the merits of both models and supported interpretation 2. Our results shed new light on this question. In our model, the power-law component can be attributed to thermal Comptonization, taking place in the hot plasma below the radiative shock. However, a natural consequence of our model is that the temperature falls off on the polar cap with the distance from the accretion column, which would be consistent with the double blackbody fit discussed by Enoto et al. (2011). It seems that at present the observational data and the model predictions are not good enough to distinguish between these alternatives or their combination.

It is an interesting question why these AXPs/SGRs have hard X-ray tails up to 100 keV or more (Enoto et al. 2010), while in general the accreting binary neutron stars show an exponential cutoff at energies of 10–20 keV. The major difference between the two classes is that AXPs/SGRs have luminosities that are about two or three orders of magnitude lower than typical X-ray binaries like Her X-1. To our knowledge, the only known accreting X-ray pulsars in the binary system, which come close to 4U 0142+61 with respect to luminosity and magnetic field strength, are 4U 0352+309 (X Persei) and 4U 2206+54, which we discuss below. Although there are a few additional sources that could be considered as members of this class of objects (e.g., 4U 1145−61 and 4U 1258−61), the scarce available hard X-ray data on them do not allow us to investigate the nature (and significance in the case of 4U 1258−61) of their hard X-ray emission at the level of detail that is possible for X Persei and 4U 2206+54.

4U 0352+309 (X-Per): according to BeppoSAX observations of di Salvo et al. (1998), this source has a luminosity of 2.4 × 10³⁴ erg s⁻¹, with half of the luminosity in the 0.1–10 keV band and the other half in the 10–100 keV band. A crude estimate derived from the spin-up rate gives a magnetic dipole field for 4U 0352+309 of B ∼ 2.5 × 10¹²(1 + z) G (Ghosh & Lamb 1979). The detection of a cyclotron line at ∼29 keV, implying a magnetic field strength of B ∼ 2.5 × 10¹²(1 + z) G, has been reported by Coburn et al. (2001). Thus, the magnetic field strength of 4U 0352+309 (X Persei) is comparable with that expected for fallback-disk accretors. The spectrum of the source exhibits a soft thermal component below 20 keV and a distinct hard tail (di Salvo et al. 1998). The photon number spectral index of the hard component between 30 and 70 keV is Γ_h ∼ 1.35, while for the hard tail of 4U 0142+61 it is Γ_h ∼ 0.93 (den Hartog et al. 2008; see also Trümper et al. 2010) or Γ_h ∼ 0.9 (Enoto et al. 2011). Thus, this source has very similar spectral characteristics as 4U 0142+61. Recently, Lutovinov et al. (2012) described the activity of the source during 2001–2011 and reported several consecutive outbursts during which the luminosities increased up to 1.2 × 10³⁵ erg s⁻¹ and a spectrum showing a distinct hard X-ray tail up to energies of ∼160 keV. The cyclotron line at 29 keV is confirmed.

We suggest that the power-law tail of 4U 0352+309 is due to the same mechanism operating in 4U 0142+61, namely, BMC.⁵

4U 2206+54: this is a pulsar showing very slow pulsations (P = 5560 s; Reig et al. 2009) accreting matter from the stellar wind of an O-type companion. A detailed analysis of the magnetorotational evolution requires a magnetic field of the neutron star B > 5 × 10¹³ G (Ikhsanov & Beskrovnaya 2011). This is in conflict with a possible cyclotron resonance scattering feature at ∼30 keV claimed by Torrejon et al. (2004), Masetti et al. (2004), Blay et al. (2005), and Wang (2009), which would indicate a magnetic field strength of ∼3.3 × 10¹² G, but which could not be detected by other observations (Reig et al. 2009; Wang 2010). At any rate, the dipole magnetic field must be large, although it will be smaller than the magnetic fields usually assigned to magnetars (Reig et al. 2012). The source shows intensity variations of a factor of five in the soft band (2–10 keV), ranging from 1.5 to 8 × 10⁻³ photons cm⁻² s⁻¹. The average luminosity in the 2–10 keV band is 1.5 × 10³⁵ erg s⁻¹ for a distance of 2.6 kpc (e.g., Reig et al. 2012). Thus, the luminosity and the magnetic dipole strength are comparable with those of AXPs/SGRs.

But there are also other similarities. (1) The soft spectrum can be modeled by a blackbody and a power-law component. The blackbody radius is ∼2.6 km, which is consistent with the size of a hot polar cap of the neutron star (Reig et al. 2012). (2) At energies >20 keV a power-law spectrum is observed that has been measured up to energies of more than 100 keV whose steepness depends on the intensity level. At low luminosity levels (≳ 10³⁴ erg s⁻¹), the source shows a hard power-law tail at energies >20 keV ranging up to more than 100 keV with a slope Γ_h ∼ 2 (Torrejon et al. 2004; Masetti et al. 2004; Reig et al. 2009). At higher luminosity levels (∼2.7 × 10³⁵ erg s⁻¹) the hard spectrum becomes significantly steeper with Γ_h ∼ 3 (Reig et al. 2012). Actually, its overall shape approaches that of high-luminosity pulsars like Her X-1, which shows a spectral cutoff at ∼20 keV.

Concluding this discussion, we note that the similarity in the spectral behavior of 4U 0352+309 and 4U 2206+54 compared with AXPs/SGRs on the one hand and the difference from the majority of the high-luminosity binary X-ray sources on the other support our view (see Section 3.2) that in the luminosity range around 10³⁵ erg s⁻¹ the transverse optical depth of the accretion column (τ ∼ 5) is favorable for the production and escape of BMC photons.

• 1.  
  We have shown in a semi-quantitative way that bulk-motion/thermal Comptonization in an accretion column, formed by a dipole magnetic field of strength B ∼ 10¹³ G, describes well not only the soft and hard X-ray spectra but also the phase-dependent energy spectra of AXP 4U 0142+61.
• 2.  
  Furthermore, we have shown in a quantitative way that the model explains naturally the observed pulse profiles in all X-ray bands. The energy-dependent pulse profiles and their constancy over long periods of time constitute significant observational constraints for the proposed models. Our model explains them in a simple and natural way by the formation of two "beams," one perpendicular to the accretion column (main pulse) and another parallel to it (secondary pulse).
• 3.  
  Our beaming model takes the gravitational bending and a noncircular shape of the accretion column into account, allowing us to determine both the angle a between the spin axis and the magnetic dipole axis and the inclination i of the spin axis with respect to the line of sight. Assuming that the radiative shock is located very close to the stellar surface (H = 0), we find a good fit with i ∼ 40° and α ∼ 20°. However, the ratio of polar to fan beam flux is too small, as only 30% of the impinging fan beam flux is converted into the scattered and thermal albedo. For a height of the radiative shock of H = 2 km, we find i ∼ 60°, α ∼ 30°, and a more reasonable value for the albedo, which is ∼75%. We conclude that our fits provide evidence for a height H of the order of at most 3.5 km.
• 4.  
  The Suzaku observations (Enoto et al. 2010) have provided a homogeneous set of spectra of two SGRs and five AXPs, showing that the fluxes of the soft F_s and the hard F_h spectral components are roughly equal, with variations of a factor of three. In view of the similarity of their broadband spectra, it is no surprise that accretion models with BB, TC, and BMC components can fit their broadband (0.5–200 keV) spectra (A. Zezas et al., in preparation). The observed ratio of these observed fluxes F_h/F_s will depend on the range of polar angles θ swept by the line of sight during one revolution. In addition, one expects in our model that the hard luminosity produced by BMC correlates with the accretion rate, and for disk accretion the latter is expected to correlate with the period derivative $\dot{P}$. Thus, one expects a correlation of F_h/F_s with $\dot{P}$. This is in agreement with the Suzaku data, which show a correlation of F_h/F_s with $\dot{P}$.
• 5.  
  We suggest that the formation of hard X-ray spectral tails in accretion columns depends crucially on the optical depth of the accretion column, which is expected to correlate with X-ray luminosity. For high-luminosity sources, the transverse optical depth τ_t of the accretion column is high, and the high-energy photons produced by BMC are thermalized before they escape. In AXPs and SGRs, having luminosities a factor of ∼100 less, the optical depth across the column is large enough to ensure efficient upscattering of photons to high energies (E ≳ 100 keV) and low enough to allow their escape before thermalization. We suggest that this requires τ_t ∼ 5.
• 6.  
  The fact that the accreting X-ray binary 4U 0352+309 (X Per), having a luminosity about a factor of 10 smaller than 4U 0142+61, shows a luminous hard tail strongly supports our model. In this case, the optical depth may be expected to be even smaller, which would explain the relatively low cutoff of the hard component (∼65 keV). The low-luminosity accreting neutron star 4U 2203+54 has a hard tail as well, showing a significant steepening of its hard spectrum at transitions from low to high luminosities, which is in line with the conclusion of point 5.
• 7.  
  It is well known that the polar caps of AXPs/SGRs show significantly higher photospheric temperatures compared with other relatively young neutron stars, like radio pulsars or X-ray dim neutron stars (XDINs) (e.g., Aguilera et al. 2008). Evidently, our model provides a natural explanation for this observational fact, since the strong illumination of the photosphere by the fan beam acts as an additional heat source. The efficiency of this process is increased by gravitational bending. The heat influx into the photosphere will depend primarily on the power of the fan beam and its width, as well as on the gravity of the neutron star.
• 8.  
  We propose that the spectral bump of the secondary pulse (which corresponds to the polar beam) is caused by cyclotron resonance reflection of the fan beam by the magnetized photosphere. This leads to magnetic field strength of B ∼ 5 × 10¹²(1 + z) G, which falls in the range of dipole magnetic fields (10¹²–10¹³ G) attributed to AXPs/SGRs in the framework of the fallback-disk scenario. Support for the above proposal is the steepness of the polar beam spectrum in the 10–60 keV range.
• 9.  
  We stress that the few giant bursts observed from SGRs cannot be explained by any accretion process. The same is true for other bursts with large super-Eddington luminosities. They are most likely caused by processes taking place at the stellar surface, e.g., conversion of super-strong magnetic fields (B ∼ 10¹⁵ G) into radiation triggered by crustal shifts, as discussed in the classical magnetar literature (e.g., Thompson & Duncan 1995). Our analysis suggests that these events do not take place in the dipole field, but in localized multipole fields. This situation is qualitatively similar to that of the Sun, which shows flare activities in sunspot fields, which are larger than the solar dipole field by at least two orders of magnitude.
• 10.  
  Finally, we note that the kind of analysis presented in this paper—applied to other AXPs/SGRs observed with sufficiently large photon statistics—can be used to study in detail the radiative properties of accretion columns and to narrow the gap between observations and theoretical models in this important field of astrophysics.

This research has been supported in part by EU Marie Curie project no. 39965, EU REGPOT project number 206469, and EU FPG Marie Curie Transfer of Knowledge Project ASTRONS, MKTD-CT-2006-042722. Ü.E. acknowledges research support from TÜBİTAK (The Scientific and Technical Research Council of Turkey) through grant 110T243 and support from the Sabancı University Astrophysics and Space Forum. We also thank Xiao-Ling Zhang (MPE) for her help with a search of the ISGRI data bank.

Facilities: INTEGRAL - International Gamma-Ray Astrophysics Laboratory satellite, XMM - Newton X-Ray Multimirror Mission satellite