A UNIFIED TIMING AND SPECTRAL MODEL FOR THE ANOMALOUS X-RAY PULSARS XTE J1810−197 AND CXOU J164710.2−455216

In recent years, an increasing number of high-resolution spectral and timing observations of isolated neutron stars has become available. Many of these observations concern two peculiar classes of high-energy pulsars, the anomalous X-ray pulsars (AXPs: nine objects plus one candidate) and the soft gamma-ray repeaters (SGRs: six objects).⁴ Historically these two classes of sources were regarded as distinct. While SGRs were first discovered in late 1978–early 1979, when SGR 1806−20 and SGR 0526−66 exhibited a bright burst of soft γ-rays (Mazets et al. 1979; Laros et al. 1986), AXPs were observed for the first time in 1981, when Fahlman & Gregory (1981) discovered pulsations in the EINSTEIN source 1E 2259+586. It was, however, not until the mid 1990s that AXPs were recognized as a class of "anomalous" pulsars because of their luminosity substantially exceeding rotational energy losses (Mereghetti & Stella 1995).

Although SGRs were mainly known as emitters of short, energetic bursts, they are also persistent X-ray sources with properties quite similar to those of AXPs (see, e.g., Woods & Thompson 2006; Mereghetti 2008, for reviews). They all are slow X-ray pulsars, with spin periods in a very narrow range (P ∼ 2–12 s), relatively large spin-down rates ($\dot{P} \sim 10^{-13}$–10⁻¹¹ s s⁻¹), spin-down ages of 10³–10⁴ yr, and stronger magnetic fields compared to those of rotation or accretion powered pulsars (B ∼ 10¹⁴–$10^{15} \ {\rm G} > B_{\rm{QED}} \simeq 4.4 \times 10^{13} \ {\rm G}$). AXPs and SGRs have persistent X-ray luminosities L_X ∼ 10³⁴–10³⁶ erg s⁻¹. Their spectra in the 0.1–10 keV band are relatively soft and can be empirically fitted with a two-component model, an absorbed blackbody (kT ∼ 0.2–0.6 keV) plus a power law (Γ ∼ 2–4). International Gamma-Ray Astrophysics Laboratory observations revealed the presence of sizeable emission up to ∼200 keV, which accounts for up to 50% of the total flux. Hard X-ray spectra are well represented by a power law, which dominates above ∼20 keV in AXPs.

The large high-energy output cannot be explained in terms of rotational energy losses, as in conventional models for radio pulsars, while the lack of stellar companions argues against accretion. The powering mechanism of AXPs and SGRs, instead, is believed to reside in the neutron star ultra-strong magnetic field (magnetar; Duncan & Thompson 1992; Thompson & Duncan 1993). The magnetar scenario appears capable of explaining the properties of both the bursts (Thompson & Duncan 1995) and the persistent emission (the twisted magnetosphere model; Thompson et al. 2002; Zane et al. 2009, and references therein; see Section 3 for details), although no definite model for the hard tails has been put forward as yet.

The persistent emission of SGRs and AXPs is now known to be variable. AXPs, in particular, display different types of X-ray flux variability: from slow, moderate flux changes on timescales of months/years, to intense outbursts with short rise times (∼1 day) lasting ∼1 yr. Some AXPs were found to undergo intense and dramatic SGR-like burst activity on sub-second timescales (XTE J1810−197, 4U 0142+614, 1E 1048.1−5973, CXOU J164710.2−455216, and 1E 2259+586). The discovery of bursts from AXPs is regarded as further evidence in favor of a common nature of AXPs and SGRs.

The first case of AXP flux variability was observed in 2002, when 1E 2259+586 showed an increase in the persistent flux by a factor ∼10 with respect to the quiescent level, followed by the emission of ∼80 short bursts with luminosity L_X ∼ 10³⁶–10³⁸ erg s⁻¹ (Kaspi et al. 2003). In early 2003, the 5.54 s AXP XTE J1810−197 was discovered at a luminosity ∼100 greater than its quiescent value (10³³ erg s⁻¹; Ibrahim et al. 2004). Analysis of archival data revealed that the outburst started between 2002 November and 2003 January.

On 2006 September 21, an outburst was observed from the AXP CXOU J164710.2−455216 (P = 10.61 s). The flux level was ∼300 times higher than that measured only five days earlier by XMM-Newton (Muno et al. 2006b; Campana & Israel 2006; Israel & Campana 2006). This, much as in the case of XTE J1810−197, indicates that some AXPs are transient sources (dubbed transient AXPs, TAXPs) and may become visible only when they enter an active state. Recently, other AXPs and SGRs showed a series of short bursts of soft γ-rays which was detected by different satellites (Mereghetti et al. 2009).

In this paper, we present a comprehensive study of the pulse profile and spectral evolution of the TAXPs XTE J1810−197 and CXOU J164710.2−455216 throughout their outbursts of 2002 November and 2006 September, respectively. By confronting timing data with synthetic light curves obtained from the twisted magnetosphere model (Nobili et al. 2008), we were able to estimate how the physical parameters of the source (surface temperature and emitting area, electron energy, twist angle) evolve in time. The fits of the pulse profiles also allowed us to infer the geometry of the two systems, i.e., the angles between the magnetic and rotational axes and the line of sight (LOS). Spectral models, obtained with the parameter values derived from the timing analysis, provide acceptable fits to XMM-Newton data.

2.1. XTE J1810−107

2.2. CXOU J164710.2−455216

It is now widely accepted that AXPs and SGRs are magnetars, and that their burst/outburst activity, together with the persistent emission, is powered by their huge magnetic field. In particular, the soft X-ray spectrum (∼1–10 keV) is believed to originate in a "twisted" magnetosphere (∇ × B ≠ 0; Thompson et al. 2002), where the currents needed to support the field provide a large enough optical depth to RCS of thermal photons emitted by the star surface. Since charges are expected to flow along the closed field lines at relativistic velocities, photons gain energy in the (resonant) scatterings and ultimately fill a hard tail.

Most studies on spectral formation in a twisted magnetosphere (Lyutikov & Gavriil 2006; Fernandez & Thompson 2007; Nobili et al. 2008) are based on the axially symmetric, force-free solution for a twisted dipolar field presented by Thompson et al. (2002). This corresponds to a sequence of magnetostatic equilibria which, once the polar strength of the magnetic field B_p is fixed, depends only on a single parameter: the radial index of the magnetic field p (B ∝ r^−p−2, 0 ⩽ p ⩽ 1) or, equivalently, the twist angle

$$\begin{equation} \Delta \phi _{N-S} = \lim _{\theta _{0} \rightarrow 0}2\int _{\theta _{0}}^{\pi /2} \frac{B_{\phi }}{B_{\theta }}\frac{d\theta }{\sin \theta }\,, \end{equation} \tag{ 1 }$$

where B_r, B_θ, and B_ϕ are the spherical components of the field, which depend only on r and θ because of axial symmetry. Knowledge of B fixes the current density j = (c/4π)∇ × B, and, if the particle velocity is known, also the electron density in the magnetosphere

$$\begin{equation} n_{e}(r,\theta) = \frac{p+1}{4\pi e}\Bigl (\frac{B_{\phi }}{B_{\theta }}\Bigr) \frac{B}{r |\langle \beta \rangle |}, \end{equation} \tag{ 2 }$$

where e is the electron charge and 〈β〉 is the average charge velocity (in units of c). The charge density of the space charge-limited flow of ions and electrons moving along the closed field lines is orders of magnitude larger than the Goldreich–Julian density, $n_{\rm{GJ}}$, associated with the charge flow along the open field lines in radio pulsars.

In our investigation, we make use of the spectral models presented by Nobili et al. (2008, NTZ in the following), who studied radiative transfer in a globally twisted magnetosphere by means of a three-dimensional Monte Carlo code. Each model is characterized by the magnetospheric twist Δϕ_N−S, the electron (constant) bulk velocity β, and the seed photon temperature kT. The polar field was fixed at B_p = 10¹⁴ G. In the applications presented by NTZ, it was assumed that the star surface emits unpolarized, blackbody radiation and is at a uniform temperature. Concerning the present investigation, the most critical assumption is that of a globally twisted magnetosphere, as discussed in more detail in Section 5. Taking a constant value for the electron bulk velocity is certainly an oversimplification and reflects the lack of a detailed model for the magnetospheric currents. In a realistic case, one would expect that β is a function of position. However, resonant scattering is possible only where the bulk velocity is mildly relativistic. If along a flux tube there are large variations of the Lorentz factor, only the region where β ≈ 0.5 will contribute to scattering. Moreover, preliminary calculations of the dynamics of charged particles in a twisted force-free magnetosphere performed accounting for both electrostatic acceleration and Compton drag indicate that β is indeed fairly constant along the central part of a flux tube (A. M. Beloborodov 2010, private communication). The assumption of unpolarized thermal radiation is not cogent either, since we are not interested in the polarization of the escaping radiation and the emergent spectrum is quite insensitive to the polarization fraction of the seed photons (see, e.g., Figure 4 of NTZ).

The code works by dividing the stellar surface into N_Θ × N_Φ zones of equal area by means of a (cos Θ, Φ) grid, where Θ is the magnetic colatitude and Φ is the longitude. After a few scatterings, photons escape from the neutron star magnetosphere and are collected on a spherical surface (the "sky") which is divided into $N_{\Theta _s} \times N_{\Phi _s}$ patches, similarly to what is done for the star surface. The key point is that the evolution of seed photons from each patch is followed separately. This allows us to treat an arbitrary surface temperature distribution without the need to perform new Monte Carlo runs, by simply combining together models from runs with different temperatures at the post-production level (the geometry is shown in Figure 1).

Zoom In Zoom Out Reset image size

Monte Carlo models are computed (and stored) for the simplest geometrical case, in which the spin and the magnetic axes are aligned. As discussed in NTZ, the most general situation in which the spin and magnetic axes are at an arbitrary angle ξ can be treated at the post-production level. If χ is the inclination of the LOS with respect to the star spin axis and α is the rotational phase angle, the coordinates of the points where the LOS intersects the sky can be found in terms of ξ, χ, and α. The pulse profile in any given energy band is then obtained by integrating over the selected range the energy-dependent counts at these positions as the star rotates (see again NTZ for details). In order to compare model light curves with observations, integration over energy is performed by accounting for both interstellar absorption and the detector response function. Actually, the interstellar absorption cross section σ and the response function A depend on the photon energy at infinity $\bar{E}=E\sqrt{1-R_S/R_{NS}}$, where E is the energy in the star frame (which is used in the Monte Carlo calculation) and R_S is the Schwarzschild radius (we assume a Schwarzschild spacetime and take R_NS = 10 km and M_NS = 1.44 M_☉). Our model pulse profile in the $[\bar{E}_1,\, \bar{E}_2]$ energy band is then proportional to

$$\begin{equation} \int _{\bar{E}_1}^{\bar{E}_2}d\bar{E}\, \exp {[-N_H\sigma (\bar{E})]}A(\bar{E}) N(\alpha, E), \end{equation} \tag{ 3 }$$

where N_H is the hydrogen column density and N(α, E) is the phase- and energy-dependent count rate. In the applications below we used the Morrison & McCammon (1983) model for interstellar absorption and, since we deal with XMM-Newton observations, we adopted the EPIC-pn response function. We remark that the Monte Carlo spectral calculation is carried out assuming a flat spacetime (i.e., photons propagate along straight lines), so that, apart from the gravitational redshift, no allowance is made for general-relativistic effects (see Zane & Turolla 2006, for a more detailed discussion). In particular, no constraints on the star mass and radius can be derived in the present case from the comparison of model and observed pulse profiles (see, e.g., Leahy et al. 2008, 2009).

Finally, phase-averaged spectra are computed by summing over all phases the energy-dependent counts. Note that 0 ⩽ ξ ⩽ π/2, while χ is in the range [0, π] because of the asymmetry between the north and south magnetic poles introduced by the current flow.

Our first step in the study of the two TAXPs XTE J1810−197 and CXOU J164710.2−455216 was to reproduce the pulse profiles (and their time evolution) within the RCS model discussed in Section 3. The fit to the observed pulse profiles in different energy bands (total: 0.5 keV ⩽ E ⩽ 10 keV, soft: 0.5 keV ⩽ E ⩽ 2 keV, hard: 2 keV ⩽ E ⩽ 10 keV) provides an estimate of the source parameters, including the two geometrical angles ξ and χ. While the twist angle, electron velocity, and surface temperature may vary in the different observations (although they must be the same in the different energy bands for a given observation), the fits have to produce values of ξ and χ which are at all epochs compatible with one another (to within the errors) in order to be satisfactory. We then computed the phase-averaged spectra for the two sources at the various epochs for the same sets of parameters and compared them with the observed ones. There are several reasons which led us to choose such an approach. The main one is that, as discussed in NTZ (see also Zane et al. 2009), spectral fitting alone is unable to constrain the two geometrical angles. Moreover, light curve fitting allows for a better control in the case in which the surface thermal map is complex and changes in time (see below).

For the present investigation, a model archive was generated beforehand. Each model was computed by evolving N_patch = 225,000 photons for N_Θ × N_Φ = 8 × 4 = 32 surface patches (N_tot = 7, 200,000 photons). The parameter grids are 0.1 ⩽ kT ⩽ 0.9 keV (step 0.05 keV), 0.1 ⩽ β ⩽ 0.9 (step 0.1), and 0.2 rad ⩽ Δϕ_N−S ⩽  1.2 rad (step 0.1 rad). Photons are collected on a $N_{\Theta _{s}} \times N_{\Phi _{s}} = 10 \times 10 = 100$ angular grid on the sky, and in N_E = 50 energy bins, equally spaced in log E in the range 0.1–100 keV.

The analysis proceeds as follows. We first used the principal component analysis (PCA) to explore the properties of the light curves as a population and to select the model within the archive that is closest to the observed one at a given epoch. This serves as the starting point for the pulse profile fitting procedure, which we performed by assuming that the whole star surface is at the same temperature. The fitting is then repeated first for the case in which the surface thermal distribution consists of a hot spot and a cooler region, and then by generating a new archive with a finer surface gridding, and applying it in the case of a surface thermal map consisting of a hot spot, a warm corona, and a cooler region (see again Figure 1). Finally, the source parameters derived from the light curve fitting are used to confront the model and observed (phase-averaged) spectra. Phase-resolved spectral analysis, although feasible in our model and potentially important, was not attempted because the decay in flux of both sources makes the counting statistics rather poor after the first one/two observations (see Bernardini et al. 2009, for more details in the case of XTE J1810−197).

4.1. PCA

The PCA is a method of multivariate statistics that allows us to reduce the number of variables X_i needed to describe a data set by introducing a new set of variables, the principal components (PCs) Z_i. The PCs are linear combinations of the original variables and are such that Z₁ displays the largest variance, Z₂ the second largest, and so on. By using the PCs it is possible to describe the data set in terms of a limited number of variables, which, however, carry most of the information contained in the original sample (see, e.g., Zane & Turolla 2006, and references therein).

Synthetic light curves were generated for 32 phases in the range [0, 2π] and for a 9 × 9 angular grid, 0° ⩽ ξ ⩽ 90° (step 10°), 0° ⩽ χ ⩽ 180° (step 20°); the archive contains a total of 136,323 models. Once the PCA was applied to the light curve set, we found that the first three PCs (Z₁, Z₂, Z₃) account for as much as ∼90% of the sample variance. This means that the entire set is satisfactorily described in terms of just three variables instead of the original 32 (see Zane & Turolla 2006, for an interpretation of Z₁, Z₂, Z₃). A graphic representation of the light curves in the archive in terms of the first three PCs is shown in Figure 2. In the same plot, we also show the PC representation of the pulse profiles of XTE J1810−197 and CXOU J164710.2−455216 at the various epochs. The points corresponding to observations fall within the volume occupied by models and this guarantees that there is a combination of the parameters for which a synthetic pulse profile reproduces the data. The PC representation is also used to find the model in the archive which is closest to a given observed light curve, by looking for the minimum of the (squared) Euclidean distance ∑³²_i=1(Z_i − Z^obs_i)² between the model and the observed pulse profile.

Zoom In Zoom Out Reset image size

4.2. XTE J1810−197

We considered eight XMM-Newton observations, covering the period 2003 September–2007 September for the TAXP prototype XTE J1810−197 (see Table 1 for the observation log). Only EPIC-pn data were used, and we refer to Bernardini et al. (2009), who analyzed the same observations, for all details on data extraction and reduction. All the EPIC-pn spectra were rebinned before fitting, to have at least 40 counts per bin and prevent oversampling the energy resolution by more than a factor of 3.

Table 1. XTE J1810−197 XMM-Newton Observations^a

Label	ObsID	Epoch (yyyy-mm-dd)	Exposure Time (s)	Total Counts	Background Counts
Sep03	0161360301	2003-09-08	5199	60136	2903
Sep04	0164560601	2004-09-18	21306	89082	1574
Mar05	0301270501	2005-03-18	24988	54279	1760
Sep05	0301270401	2005-09-20	19787	21876	1311
Mar06	0301270301	2006-03-12	15506	12296	1197
Sep06	0406800601	2006-09-24	38505	23842	2974
Mar07	0406800701	2007-03-06	37296	21903	2215
Sep07	0504650201	2007-09-16	59014	34386	4117

Note. ^aEPIC-pn

Download table as:  ASCIITypeset image

4.2.1. Pulse Profiles

We started our analysis by making the simplest assumption about the star surface thermal map, a uniform distribution at temperature T. Light curves were then computed in the total, soft, and hard energy band for all the models in the archive. Once the model closest to each observation (and in each band) was found through the PCA, we used it as the starting point for a fit performed using an IDL script based on the minimization routine mpcurvefit.pro. Our fitting function has six free parameters, because, in addition to the twist angle, the temperature, the electron velocity, and the angles χ and ξ, we have to include an initial phase to account for the indetermination in the position of the pulse peak. Since it is not possible to compute "on the fly" the pulse profile for a set of parameters different from those contained in the archive, light curves during the minimization process were obtained from those in the archive using a linear interpolation in the parameter space.

In this way, we obtained a fair agreement with the observed pulse profiles (χ² ⩽ 1.12 in five out of eight observations; see Table 4), and the values of the physical parameters (Δϕ_N−S, β, T) turn out to be the same (to within the errors) for a given epoch among the different energy bands, as they need to be. Moreover, the evolution of the twist angle and of the surface temperature follows a trend in which both quantities decrease in time as the outburst declines. This is expected if the outburst results from a sudden change in the NS magnetic structure, producing both a heating of the star surface layers and a twisting of the magnetosphere which then dies away (Thompson et al. 2002; Beloborodov 2009). However, the model is not acceptable since we found that the geometrical angles χ and ξ change significantly from one observation to another, and even for the same observation in the different energy bands (see Figure 3 where the parameter evolution is shown for the three energy bands). The analysis of the hard band was not carried out after 2006 September, because in both of the 2007 observations there are only a few photons with energy >2 keV and, as a consequence, light curves are affected by large uncertainties.

Zoom In Zoom Out Reset image size

This shortcoming is most probably due to our oversimplifying assumption about the NS thermal map. In fact, it was shown that in the post-outburst phase the surface temperature distribution of XTE J1810−197 is complex and changes in time (Perna & Gotthelf 2008; Bernardini et al. 2009, although a different emission model was assumed in these investigations). In order to check this, we tried different configurations, starting with a two-temperature map: a hot cap centered on one magnetic pole with the rest of the surface at a constant, cooler temperature. In this picture, both temperatures and the emitting areas are allowed to vary in time. While applying the one-temperature model we found that both the 2007 observations were reasonably well reproduced with a value of the (uniform) surface temperature ∼0.15 keV, comparable to the quiescent one (see also Bernardini et al. 2009). In order to check if this fit can be further refined, we started from the 2007 September observation, freezing the colder temperature at T_c = 0.15 keV, and letting the hot cap temperature T_h free to vary. Since the cap area A_h is not known a priori, nor it can be treated as a free parameter in our minimization scheme, we tried several values of A_h, corresponding to one up to eight patches of our surface grid (this means that A_h is n/32 of the star surface, with n = 1,  ...,  8). The best-estimate emitting area was then taken as the one giving the lowest reduced χ² for the fit in the different trials. We verified that in all cases the same value of the cap area produces the minimum χ² in all energy bands. Independent of the emitting area chosen, we always found for T_h a value compatible with ∼0.15 keV for both the 2007 September and March observations.

One can then conclude that, for these two epochs, the entire star is radiating at the same temperature, or if a hot cap exists, its area is smaller than ∼3% of the star surface (the size of our surface grid resolution). For these epochs, we report in Table 2 the values of the cold temperature obtained by using the single temperature scenario. We note that the temperature value is at the border of our grid of parameter values, so that, strictly speaking, it should be regarded as an upper limit on T_c. However, we verified that the χ² grows steeply when T_c increases above 0.15–0.16 keV. Although there is no guarantee that the same is true when T_c decreases, in the following we assume that T_c ∼ 0.15 keV is a satisfactory estimate for the uniform temperature at these epochs.

Table 2. XTE J1810−197 Parameters and Thermal Map (Two-temperature Model)^a

Epoch	ΔϕN−S	β	ξ(°)	χ(°)	Th (keV)	Ah(%)	Tc (keV)	Ac(%)
Sep03	0.70 ± 0.01	0.80 ± 0.01	22.7 ± 0.5	144.2 ± 0.6	0.71 ± 0.01	25.0 ± 3.1	0.30	75.0 ± 3.1
Sep04	0.67 ± 0.01	0.62 ± 0.02	20.7 ± 1.0	158.2 ± 0.2	0.55 ± 0.01	18.7 ± 3.1	0.30	81.3 ± 3.1
Mar05	0.61 ± 0.01	0.49 ± 0.01	21.6 ± 0.4	147.1 ± 0.6	0.67 ± 0.01	12.5 ± 3.1	0.25	87.5 ± 3.1
Sep05	0.47 ± 0.01	0.53 ± 0.05	23.0 ± 0.1	159.0 ± 1.1	0.42 ± 0.01	9.4 ± 3.1	0.25	90.6 ± 3.1
Mar06	0.49 ± 0.01	0.50 ± 0.11	23.5 ± 0.4	149.4 ± 3.5	0.28 ± 0.01	6.2 ± 3.1	0.15	93.8 ± 3.1
Sep06	0.43 ± 0.01	0.71 ± 0.16	21.4 ± 0.3	155.7 ± 1.7	0.28 ± 0.01	3.1 ± 3.1	0.15	96.9 ± 3.1
Mar07	0.45 ± 0.01	0.60 ± 0.01	29.8 ± 0.1	162.8 ± 0.1	...	...	0.16 ± 0.01	100.0
Sep07	0.48 ± 0.01	0.70 ± 0.08	22.4 ± 0.1	163.0 ± 1.7	...	...	0.15 ± 0.01	100.0

Note. ^aTotal energy band; parameters with no reported errors are fixed. Parameter errors are calculated by the minimization routine mpcurvefit.pro and are at 1σ. Errors on the area correspond to the smallest patch of the grid.

Download table as:  ASCIITypeset image

We then proceeded backward in time, from 2006 September till 2003 September. Again, the cooler temperature is kept fixed while several values of A_h are tried. However, to account for the possibility that T_c also varies, we repeated the calculation for T_c = 0.15,  0.20,  0.25,  0.30 keV, looking for the pair (T_c,  A_h) which gives the lowest χ². Results are summarized in Table 2. Although the fits improve with respect to the one-temperature model (see Table 4), the two geometrical angles still change from one observation to another and also across different bands at the same epoch.

In order to more accurately reproduce the star thermal map, we generated a new model archive, increasing the number of surface patches to N_Θ × N_Φ = 50 × 4 = 200. The temperature, electron velocity, and twist angle are in the range 0.15 keV ⩽ T ⩽ 0.9 keV (step 0.15 keV), 0.1 ⩽ β ⩽ 0.9 (step 0.2), and 0.4 rad ⩽ Δϕ_N−S ⩽ 1.2 rad (step 0.2 rad), respectively. We then assumed that the star surface is divided into three zones: a hot cap at temperature T_h, a concentric warm corona at T_w, and the remaining part of the neutron star surface at a cooler temperature T_c. Again, we began our analysis from the 2007 observations, fixing T_c = 0.15 keV, and searching for the value of the warm temperature T_w. Every fit was repeated for 12 values of the emitting area A_w = 0.5%,  1%,  2%,  4%,  ...,  20% the total surface. We found that the reduced χ² improves with the addition of a warm cap at T_w ∼ 0.3 keV, accounting for 0.5% of the neutron star surface (see Table 3). We stress that this value is below the resolution of our previous grid, so the two results are consistent with each other.

Table 3. XTE J1810−197 Parameters and Thermal Map (Three-temperature Model)^a

Epoch	ΔϕN−S	β	ξ(°)	χ(°)	Th (keV)	Tw (keV)	Tc (keV)	Ah(%)	Aw(%)
Sep03	0.80+0.05−0.11	0.70+0.08−0.06	27.8+4.6−3.1	145.3+4.7−2.5	0.62+0.14−0.14	0.30	0.15	8 ± 0.5	16 ± 0.5
Sep04	0.79+0.07−0.08	0.78+0.09−0.23	16.2+4.2−5.9	140.8+5.8−2.6	0.49+0.03−0.22	0.30	0.15	6 ± 0.5	14 ± 0.5
Mar05	0.62+0.03−0.03	0.51+0.07−0.09	22.2+5.4−13.1	146.9+7.9−1.8	0.49+0.01−0.04	0.30	0.15	4 ± 0.5	14 ± 0.5
Sep05	0.53+0.10−0.09	0.50+0.11−0.19	21.4+10.3−20.0	154.4+13.9−9.5	0.52−−	0.30	0.15	2 ± 0.5	10 ± 0.5
Mar06	0.46+0.08−0.04	0.73+0.200.20	18.518.0−17.2	143.0+7.9−7.5	...	0.29+0.17−0.03	0.15	...	6 ± 0.5
Sep06	0.54+0.04−0.03	0.42+0.13−0.12	22.4+12.3−20.0	150.2+14.3−9.5	...	0.30−−	0.15	...	2 ± 0.5
Mar07	0.49+0.20−0.05	0.43+0.13−0.12	30.0+12.3−20.0	153.9+19.6−16.0	...	0.29−−	0.15	...	0.5 ± 0.5
Sep07	0.47+0.07−0.04	0.50+0.09−0.15	22.7+16.4−20.0	145.8+16.4−9.5	...	0.31−−	0.15	...	0.5 ± 0.5

Note. ^aTotal energy band; parameters with no reported errors are fixed. Errors are computed from the χ² curve (see the text for details) and are at 1σ. No errors are reported when they could not be calculated (flat χ² curves) and errors on the area have the same meaning as in Table 2.

Download table as:  ASCIITypeset image

We then considered the two 2006 observations; in the two-temperature model based on the previous archive, these were reasonably reproduced with T_c = 0.15 keV and T_h ∼ 0.3 keV (note that T_h for the two-temperature model corresponds to T_w in the present case). For these two observations we repeated the fit, fixing T_c at 0.15 keV while leaving T_w free to vary. The size of the emitting area was estimated by following the same procedure discussed above. We found an almost constant value, T_w ∼ 0.3 keV, between 2006 March and 2007 September, while the emitting area decreases in time. Also, we found no need for a further component at T_h at these epochs. On the other hand, results for the two-temperature case (see Table 2) show the presence of a component with temperature higher than 0.3 keV, in the period between 2003 September and 2005 September (while the cooler one varies between 0.25 and 0.30 keV). It is tempting to associate this with a transient hot cap that appears only in the first period after the outburst, superimposed on the other, longer-lived emitting zones.

To test this possibility, we re-fitted the first four observations by fixing the coldest temperature at T_c = 0.15 keV, the warmer one at T_w = 0.3 keV, and leaving only the hotter temperature free to vary. For each observation the pulse profile fits were computed for every combination of A_h and A_w chosen among the 12 values in the range 0.5%–20% introduced before, and looking for the minimum of the reduced χ². Results of the light curve fitting at different epochs are listed in Table 3 and shown in Figure 4, while the reduced χ² for the three thermal distributions is reported in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Reduced χ² for XTE J1810−197^a

Epoch	1T $\chi ^2_{\rm{red}}$	2T $\chi ^2_{\rm{red}}$	3T $\chi ^2_{\rm{red}}$	XSPEC $\chi ^2_{\rm{red}}$	T (keV)
Sep03	1.72	1.58	0.12	1.22	...
Sep04	0.66	0.42	0.36	1.93	...
Mar05	1.02	0.98	0.79	1.50	...
Sep05	1.06	0.40	0.39	1.52	0.53+0.07−0.06
Mar06	2.94	1.70	1.25	1.34	...
Sep06	0.94	0.38	0.35	1.36	0.31+0.03−0.01
Mar07	2.88	2.88	2.37	1.08	0.29+0.04−0.02
Sep07	1.12	1.12	0.96	1.29	0.31+0.01−0.01

Note. ^aFirst three columns: reduced χ² obtained from the light curves fitting for the total energy band (results for the other two bands are similar). Last two columns: reduced χ² obtained from the spectral fitting in XSPEC, and corresponding temperatures. The temperature was left free to vary only at those epochs and for those components for which the light curve analysis did not produce a unique value. Errors for the temperature are at 1σ.

Download table as:  ASCIITypeset image

A worry may arise whether the best-fitting values obtained from the minimization routine indeed correspond to absolute minima of the reduced χ². In order to check this, and visually inspect the shape of the χ² curve close to the solution, we computed and plotted the reduced χ² leaving, in turn, only one parameter free and freezing the remaining five at their best-fit values. This also allowed us to obtain a more reliable estimate of the parameter errors which were computed by looking, as usual, for the parameter change which corresponds to a 1σ confidence level (and reported in Table 3).

We found that all values obtained with the mpcurvefit.pro routine indeed correspond to minima of the reduced χ² curve, with the exception of the temperature(s), for which there are observations (or energy bands) with very flat χ² curves (see Figure 5). In particular, for the 2005 September observation the curve obtained varying T_h is flat in all the three energy bands. Also, the curves relative to T_w for the 2006 September, 2007 March, and 2007 September observations have the same problem. This can be understood by noting that in all these observations the size of the hot/warm region accounts for only <2% of the total neutron star surface: temperature changes in such a small emitting area can hardly influence the fit. In addition, for the 2005 March and 2006 March observations the reduced χ² curve relative to one of the temperatures is flat, but this occurs for only one of the three energy bands. The first case concerns the hot temperature and the soft band, the second the warm temperature and the hard band. As we discussed above, when the hot (warm) area shrinks it has little effect on the pulse profile; this shows up first in the energy band in which its emission contributes less, i.e., the soft (hard) band.

Zoom In Zoom Out Reset image size

Given these findings, we concluded that light curve analysis by itself is unable to yield a unique temperature value for the 2005 September, the 2006 September, and both of the 2007 observations. On the other hand, spectral analysis is more sensitive to temperature variations, so that it is also possible to infer a temperature value in these cases. As will be discussed in the next section, by combining the two techniques we can remove most of the uncertainties and validate the three temperature model presented so far (see Section 4.2.2 for details).

There are several physical implications than can be drawn from our model. The TAXP is seen at an angle χ = 148^°+7₋₉ with respect to the spin axis. The misalignment between the spin axis and the magnetic axis is ξ = 23^°+15₋₁₁. These values of the two angles, and the corresponding errors, are calculated from the weighted average in the three energy bands. To get a quantitative confirmation that χ and ξ do not change in time, we fitted a constant through the values of each angle as derived from the light curve fitting at the different epochs and found that the null hypothesis probability is <1%. We note that, formally, the misalignment between the spin and the magnetic axis is compatible with being zero at the 3σ level. Low values of ξ produce, however, models with pulsed fractions quite smaller than the observed ones and, despite ξ ∼ 0 might be still statistically acceptable, we regard this possibility as unlikely because the amplitude is the main feature which characterizes the pulse, as the PCA shows (see Section 4.1, the first PC, Z₁, is, in fact, directly related to the amplitude).

A scenario emerges in which, before the outburst, the NS surface radiates uniformly at a temperature T_c ∼ 0.15 keV. Soon after the burst the thermal map of XTE J1810−197 substantially changes. The region around the magnetic north pole is heated, reaches a temperature of ∼0.7 keV, and covers an area ∼8% of the total star surface. This hot spot is surrounded by a warmer corona at ∼0.3 keV, that covers a further ∼16% of the surface. During the subsequent evolution, the hot cap decreases in size and temperature until the 2006 March observation, when it becomes too small and cold to be distinguished from the surrounding warm corona. The warm region remains almost constant until 2005 September, then decreases in size, and becomes a cap in 2006 March, following the hot spot disappearance. In 2007 September (our last observation for XTE J1810−197), the warm cap is still visible, even if its area is down to only ∼0.5% of the total. The twist angle is highest at the beginning of the outburst (2003 September) and then steadily decreases until it reaches a more or less constant value around 2005 September. The electron velocity does not show large variations in time and stays about constant at β ∼ 0.5.

Synthetic and observed light curves (in the total band) are shown in Figure 6, together with the fit residuals. We note that the residuals exhibit a well-defined, oscillatory pattern at all epochs. In our scenario, this can possibly be associated with a more complicated thermal map, of which our 3T model is a first-order approximation (e.g., non-circular shape of the hotter regions, off-centering of the hot and warm areas). However, as discussed in more detail in Section 5, no further refinement of the surface thermal map will be attempted here. Since XTE J1810−197 pulse profiles are fairly sinusoidal, we can compute the pulsed fraction and its evolution in time at different energies. The comparison of model results with data is shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Spectra

In order to verify if the thermal map inferred from the pulse profile fits is reasonable, and in order to remove the uncertainties in the value of the temperature at certain epochs (see Section 4.2.1), we examined the source spectra. The goal is to check if the parameters derived from our light curve analysis (twist angle, bulk velocity, size, and temperature of the three emitting areas, and the angles χ, ξ) can also reproduce the spectral evolution of XTE J1810−197 during the outburst decay. To this end, we used the ntzang model that was implemented in XSPEC by Nobili et al. (2008, the model is not currently available in the public library, but it can be obtained from the authors upon request). The ntzang XSPEC model has the same free parameters as those used in our fits of the pulse profiles. In addition, it contains the normalization and the column density. We have a caveat that, since this XSPEC model was created by assuming that the entire star surface emits at uniform temperature, strictly speaking it is not directly suited to the present case. As an approximation, we fitted the spectra by adding together three (absorbed) ntzang models, each associated with one of the three thermal components, at temperatures T_h, T_w, and T_c, respectively. At each epoch the fit was performed by freezing Δϕ_N−S, β, T, χ, and ξ at the values derived from the fit of the light curve in the total energy band (see Section 4.2.1), while the three model normalizations (which are related to the emitting areas) were left free to vary. We also required that the column density, N_H, be the same for all three spectral components and for all epochs. Since for the 2005 September, the 2006 September, and both of the 2007 observations, the light curve analysis did not return a unique value for the hotter temperature, we also left this parameter free to vary in these four observations. In all these cases, we found that the fit converges to a value of the temperature close to the best-fitting value obtained from the light curve analysis (see Table 4). Moreover, the reduced χ² significantly worsens by varying the temperature, meaning that the spectra are much more sensitive to the presence of these components. Results are shown in Figure 8, while the reduced χ² for the fits at the various epochs are reported in Table 4. The value of the column density is found to be N_H = (7.73 ± 0.50) × 10²¹ cm^-2, compatible at the 1.5σ level with the one obtained by Bernardini et al. (2009) with the 3BB model, N_H = (6.3 ± 0.5) × 10²¹ cm^-2. We remark that, in assessing the goodness of the fits, only the normalizations of the three components (plus N_H) are free to vary; all the other model parameters are frozen at the best values obtained from the pulse profile analysis. Under these conditions, we regard the agreement of our model with observed spectra as quite satisfactory. We note that the presence of systematic residuals at high energies (above 7–8 keV) may be hinted in the fits of the three earlier observations (see Figure 8). As discussed by Bernardini et al. (2009) they may be related to a harder spectral component which is however only marginally significant (3.2σ confidence level) and quite unconstrained. Given that the high-energy residuals are comparable to (or smaller than) those of the 3BB model used by Bernardini et al. (2009), we conclude that a hard tail is also not significant in our modeling and we did not attempt to include it in our fits.

Zoom In Zoom Out Reset image size

We checked how the reduced χ² for the spectral fit changes when the (frozen) parameters are varied within ∼2σ from their best-fit value (as from the pulse fitting). This has been done by changing one parameter at a time. We found that the χ² indeed increases quite smoothly in response to the change of each parameter, with the exception of χ and ξ. This is not surprising, since we knew already that the spectrum is not very sensitive to the geometry. We also tried a fit leaving all the parameters free, apart from the two geometrical angles which were held fixed at their best-fit values. The fit returns parameter values which are the same, within the errors, as those derived from the pulse fitting and comparable values of the reduced χ², implying that the solution we presented is indeed a global χ² minimum. The same procedure and the same conclusions also hold in the case of CXOU 164710.2−455216 (see Section 4.3.2).

4.3. CXOU J164710.2−455216

Having verified that our model can provide a reasonable interpretation for the post-outburst timing and spectral evolution of TAXP prototype XTE J1810−197, we applied it to CXOU J164710.2−455216, the other transient AXP for which a large enough number of XMM-Newton observations covering the outburst decay are available (see Table 5 for details). During 2006 September, the pn and MOS cameras were set in full window imaging mode with a thick filter (time resolution = 73.3 × 10⁻³ s and 2.6 s for the pn and MOS, respectively), while all other observations were in a large and small window imaging mode with a medium filter (time resolution = 4.76 × 10⁻² s and 0.3 s for the pn and MOS, respectively). To extract more than 90% of the source counts, we accumulated a one-dimensional image and fitted the one-dimensional photon distribution with a Gaussian. Then, we extracted the source photons from a circular region of radius 40''(smaller than the canonical 55'', corresponding to 90% of the source photons, in order to minimize the contamination from nearby sources in the Westerlund 1 cluster) centered at the Gaussian centroid. The background for the spectral analysis was obtained (within the same pn CCD where the source lies and a different CCD for the MOS) from an annular region (inner and outer radii of 45'' and 65'', respectively) centered at the best source position. In the timing analysis, the background was estimated from a circular region of the same size as that of the source. EPIC-pn spectra were processed as in the case of XTE J1810−197 (see Section 4.2).

Table 5. CXOU J164710.2−455216 XMM-Newton observations^a

Label	ObsID	Epoch (yyyy-mm-dd)	Exposure Time (s)	Total Counts	Background Counts
Sep06	0311792001	2006-09-22	26780	56934	1709
Feb07	0410580601	2007-02-17	14740	18734	1264
Aug07	0505290201	2007-08-19	16020	11710	2384
Feb08	0505290301	2008-02-15	9080	4618	1131
Aug08	0555350101	2008-08-20	26360	7357	1689
Aug09	0604380101	2009-08-24	33030	4974	1959

Note. ^aEPIC-pn

Download table as:  ASCIITypeset image

4.3.1. Pulse Profiles

The analysis of the pulse profiles of CXOU J164710.2− 455216 closely follows that presented in Section 4.2.1. In particular, we first tried a single temperature and then a two-temperature model, encountering the same problems we found for XTE J1810−197. Finally, we applied a three-zone thermal map and this provided reasonable fits for the light curves, and the angles χ and ξ were not found to vary in the same observation for the different energy bands and for different epochs. We did not attempt to fit the pulse profiles in the hard band after the 2007 February observation because of the very few counts at energies >2 keV. As for XTE J1810−197 we started the analysis from the last observation (2009 August) assuming a thermal map comprising a hot cap centered on the magnetic pole at temperature T_h, a concentric warm corona at T_w, and the rest of the neutron star at the colder temperature T_c. Every fit was repeated for 10 values of the hot cap area A_h = 0.5%,  1%,  2%, 4 %,  ...,  16% (of the total surface) and for 20 values of the warm corona area A_w = 0.5%,  1%,  2%, 4 %,  ...,  30%. Moreover, light curve fits were iterated for two values of the cold temperature T_c = 0.15, 0.30 keV and also for two values of the warm temperature T_w = 0.30, 0.45 keV. The hotter temperature was left free to vary. We found that in the last two observations, independent of the hot cap size, T_h is always ∼0.45 keV, nearly indistinguishable from the temperature of the warm corona obtained from the fit. We concluded that, at least for our present surface grid resolution, in the last two observations there are only two thermal components that contribute to the emission, the cold and warm ones, and repeated the fit leaving T_w free to vary. Results are reported in Table 6 and plotted in Figure 9, while a comparison of the reduced χ² for the three thermal distributions is given in Table 7. Errors listed in the tables have the same meaning as in the case of XTE J1810−197. Again, when the spot at T_h becomes very small its temperature cannot be determined unambiguously.

Zoom In Zoom Out Reset image size

Table 6. CXOU J164710.2−455216 Parameters and Thermal Map (Three-temperature Model)^a

Epoch	ΔϕN−S	β	ξ(°)	χ(°)	Th (keV)	Tw (keV)	Tc (keV)	Ah(%)	Aw(%)
Sep06	1.12+0.08−0.14	0.18+0.03−0.03	83.5+1.0−1.4	20.8+0.1−0.5	0.70+0.20−0.11	0.45	0.15	8 ± 0.5	22 ± 0.5
Feb07	1.07+0.05−0.10	0.19+0.02−0.03	80.0+2.7−2.1	23.2+1.7−1.7	0.64+0.16−0.06	0.45	0.15	6 ± 0.5	24 ± 0.5
Aug07	1.00+0.18−0.06	0.15+0.02−0.05	82.1+1.1−1.1	19.2+2.7−1.3	0.63+0.18−0.10	0.45	0.15	4 ± 0.5	26 ± 0.5
Feb08	0.77+0.21−0.12	0.20+0.05−0.09	85.7+2.7−4.1	20.4+7.2−0.7	0.62−−	0.45	0.15	2 ± 0.5	28 ± 0.5
Aug08	0.65+0.12−0.07	0.70+0.05−0.05	80.1+2.7−9.7	28.4+4.9−7.3	...	0.49+0.02−0.05	0.15	...	30 ± 0.5
Aug09	0.55+0.11−0.10	0.79+0.06−0.06	87.0+9.7−10.9	25.5+4.9−10.	...	0.46+0.05−0.05	0.15	...	30 ± 0.5

Note. ^aTotal energy band. Errors have the same meaning as in Table 3.

Download table as:  ASCIITypeset image

Table 7. Reduced χ² for CXOU J164710.2−455216^a

Epoch	1T χ2red	2T χ2red	3T χ2red	XSPEC χ2red	T (keV)
Sep06	1.05	0.86	0.31	1.24	...
Feb07	1.32	0.76	0.65	0.83	...
Aug07	0.97	0.91	0.44	1.01	...
Feb08	1.45	1.12	0.63	1.08	0.62+0.06−0.09
Aug08	1.45	1.23	0.79	1.23	...
Aug09	2.03	1.97	1.52	1.36	...

Note. ^aSame as in Table 4.

Download table as:  ASCIITypeset image

According to our model, CXOU J164710.2−4552116 is viewed at an angle χ = 23^°+4₋₃ with respect to its spin axis. The spin and the magnetic axes are almost orthogonal, ξ = 84^°+5₋₃. This is a quite peculiar condition, and it seems to be the only one capable of explaining the characteristic three-peaked shape of the observed light curves within the present model. As for XTE J1810−197 values and errors for both angles are calculated as the weighted average of parameters in the three energy bands. Also in this case the probability that χ and ξ are not constant in time is <1%.

Soon after the burst the thermal map of CXOU J164710.2−455216 consists of three regions at different temperatures. The hottest region, around the north magnetic pole, has a temperature T_h ∼ 0.7 keV, and its area is ∼8% of the total. This hot spot decreases in temperature and size as time elapses, until 2008 February. In 2008 August, the hot cap becomes so small in size and its temperature so close to that of the warm corona, that it is impossible to distinguish between the two regions. The warm corona has a temperature of ∼0.45 keV, which remains about constant during the three years of observations. In this case the corona area slightly increases with time, starting from ∼20% and reaching ∼30% of the NS surface. The third region has a lower temperature T_c ∼ 0.15 keV and its area remains constant at ∼70% of the total. The twist angle is ∼1.12 rad soon after the burst, and it decreases with time. There are hints that its decay is slower until 2007 August, then proceeds faster. The electron velocity is about the same at all epochs (β ∼ 0.2), apart from the last two observations in which it strongly increases. This variation may be related to the change in the pulse shape (from three-peaked to single-peaked) and also to the increase of the pulsed fraction.

A comparison between the observed and model pulse profiles is shown in Figure 10. Because of the inherent complexity and drastic time evolution of the CXOU J164710.2−455216 light curves, the agreement is not as good as for XTE J1810−197. The fact that light curve fits return χ² values not much higher than those of XTE J1810−197 (compare Tables 4 and 7) reflects the larger uncertainties in the phase-binned source counts. We also note that the errors on the geometrical angles for CXOU J164710.2−455216 are smaller than those derived for XTE J1810−197, despite the worse agreement (see Tables 3 and 6). This is most probably due to the different shapes of the pulses in the two sources. Because of the very peculiar light curve of CXOU J164710.2−455216, which can be reproduced by our model only invoking a nearly orthogonal rotator, even small departures of ξ and χ from their best-fit values results in a rapid growth of the χ². This does not occur for the rather sinusoidal pulse of XTE J1810−197 since the model can produce light curves of more or less the right shape in a wider range of angles.

Zoom In Zoom Out Reset image size

Besides being of limited use because of the complex shape of the pulse, the pulsed fraction analysis was hindered by the lower count rate, especially at low energies and was not pursued further for this source. As in XTE J1810−197, we checked that the values obtained from the minimization routine indeed correspond to minima of the reduced χ². Again we froze five of the six parameters to the value obtained with the mpcurvefit.pro minimization routine and calculated the reduced χ² around its minimum by varying the free parameter. The procedure was repeated for all parameters and all observations in the three energy bands. Again, for all parameters but the temperature, results obtained with the mpcurvefit.pro routine indeed correspond to the minima of the reduced χ² curve. There is one observation for which the χ² curve relative to the hot temperature is very flat for all the energy bands. This is the 2008 August observation, for which the size of the emitting area accounts for just 2% of the total neutron star surface. As in XTE J1810−197, we conclude that the fit is not very sensitive to the temperature variation for very small emitting areas. On the other hand, like in the previous case, it was possible to infer a value for the 2008 August hot temperature using the spectral analysis (see Section 4.3.2).

4.3.2. Spectra

The spectral analysis for CXOU J164710.2−455216 was carried out using the same approach discussed in Section 4.2.2. We fitted three ntzang components, each representative of an emitting region at different temperature, and froze all parameters apart from the three normalizations and N_H (which was forced to be the same for all the components and for all epochs). Moreover, since the light-curve analysis of the 2008 February observation failed to provide an unambiguous value for the hot temperature, this parameter was also left free to vary. Results are shown in Figure 11. Given the approach we used for the fit, the agreement is quite satisfactory (reduced χ² are listed in Table 7). Systematic residuals at low (1–2 keV) energies are however present, especially in the 2006 September, 2007 August, and 2008 August observations. N_H is found to be (2.14 ± 0.015) × 10²² cm^-2, somewhat higher than that derived by Naik et al. (2008), N_H = (1.73 ± 0.03) × 10²² cm^-2.

Zoom In Zoom Out Reset image size

The simultaneous study of the timing and spectral characteristics of the transient AXPs XTE J1810−197 and CXOU J164710.2−455216 presented in this paper shows that the post-burst evolution of two sources shares a number of similar properties. In particular, the long-term variability of the pulse profiles and spectra appears to be (semi)quantitatively consistent with a scenario in which the star surface thermal distribution and magnetospheric properties progressively change in time. Our results were derived within the twisted magnetosphere model for magnetars and support a picture in which the twist affects only a small bundle of closed field lines around one of the magnetic poles. As discussed by Beloborodov (2009), if the twist is initially confined along the magnetic axis, the returning currents hit a limited portion of the star surface (typically a polar cap), which becomes hotter. In this scenario, the post-outburst evolution is related to the twist decay, during which the bundle shrinks, and the heated region decreases both in size and temperature. We found evidence for a cooling/shrinking of the heated polar cap in both sources, together with a decrease of the twist angle. It should be noted that our magnetospheric model assumes a global twist, since no spectral calculations are currently available for a localized twist.

Within this common framework, there are nonetheless differences between the two TAXPs. For XTE J1810−197 we found that the star thermal map comprises three regions: a hot cap, a surrounding warm corona, and the rest of the surface at a colder temperature. The hot cap decreases in size and temperature until it becomes indistinguishable from the corona around 2006 March. Also, the warm corona shrinks, although its temperature stays constant at about ∼0.3 keV. It becomes a cap in 2006 March and it is still visible in our last observation (2007 September) although its size is down to 0.5% of the entire surface. The rest of the surface remains at a temperature comparable to the quiescent one (as measured by ROSAT) during the entire evolution, indicating that the outburst likely involved only a fraction of the star surface. Bernardini et al. (2009) obtained similar results using a 3BB model, although they did not attempt to locate the different emitting regions on the star surface nor to fit the pulse profiles. In their (purely spectral) analysis, the hot region is visible slightly longer (until 2006 March); the reason for the difference with respect to our results being most probably the resolution of our surface grid. Moreover, in our case the hot temperature decrease is more pronounced. The twist angle decreases from ∼0.8 rad to ∼0.5 rad during the first two years, and then it remains roughly constant.

Much as in the case of XTE J1810−197, the thermal map of CXOU J164710.2−455216 is well reproduced by three different regions. However, while the evolution of the hot cap is similar, i.e., it decreases in size and temperature until it disappears in the 2008 August observation, the behavior of the warm corona is different. Now the warm temperature remains constant at ∼0.45 keV and the area increases. Actually, the area of the "hot+warm" region is constant and covers about ∼30% of the surface, while the remaining ∼70% is at a constant cooler temperature, ∼0.15 keV. This is suggestive of a picture in which the "quiescent" state of the source is characterized by a two-temperature map, with a warm polar region superposed to the cooler surface. The outburst might have heated a portion of the warm cap, producing the hot zone which then cooled off. It is intriguing to note that the disappearance of the hot spot occurs at the same time (2008 August) at which the pulse profile dramatically changed, switching from a three-peaked to a single-peaked pattern. A quasi-sinusoidal shape of the light curve was observed when the source was in quiescence Israel et al. (2007). However, at that time the pulsed fraction was nearly 100% above 4 keV, likely indicating the presence of a small hot spot which is periodically occulted as the star rotates. This is in agreement with our finding that this TAXP is a nearly orthogonal rotator. Whether CXOU J164710.2−455216 is presently approaching quiescence is unclear. If this is the case, its quiescent state is different from that observed in 2005 and also from that of XTE J1810−197.

It is worth stressing that our claim that the temperature does not change spatially in each of the regions should not be taken literally. The assumption that the surface can be divided into three (or two) thermal regions was mainly introduced to simplify the calculations while catching the essential features of the model. A smooth temperature variation within a zone is likely to be present. However, it is difficult to reconcile the observed pulsed fraction of XTE J1810−197 in the 2006 September observation (≳10%, see Figure 7) even accounting for the temperature gradient induced by the large-scale dipolar field. This may be an indication that, as our analysis shows, there is a residual twist even in the quiescent state.

In this respect, we note that our spectral calculation is based on a rather fine subdivision of the star surface (50 × 4 patches in the final version of the archive), so we could have produced pulse profiles for arbitrary complicated thermal maps. The motivation of our choice of the thermal distribution (a hot polar cap and a warm corona superimposed to the colder surface) is threefold: (1) a model based on two thermal components, originating from a hot cap and a warm corona, was successfully applied to XTE J1810−197 by Perna & Gotthelf (2008); (2) inclusion of a third, colder component in the spectrum of the same source was shown to be statistically significant by Bernardini et al. (2009); and (3) it is consistent with theoretical predictions for a twisted magnetosphere in an AXP (Beloborodov 2009, see above). In addition, this is the simplest map for which we were able to obtain constant values, to within the errors, for the two geometrical angles χ and ξ during the entire period covered by the observations.

In their analysis of XTE J1810−197, Perna & Gotthelf (2008) assumed that the X-rays come from two concentric regions with varying temperatures and areas, each emitting a blackbody spectrum; the rest of the surface was taken to be at zero temperature. They derived the angles χ and ξ, and, although their solution is not unique, they claim that the pair χ ∼ 53°, ξ ∼ 23° is favored. While this value of ξ coincides with our estimate, the two values of the inclination of the LOS are in substantial disagreement. Also the emitting areas of the hot/warm region and their temperatures turn out to be different in the two cases. Their estimate of the hot temperature is always higher than ours and the size of the warm corona is not monotonically decreasing. We remark that quantitative differences are to be expected given the different assumed spectral models (blackbody versus RCS); moreover, because Perna & Gotthelf (2008) did not include a colder region.⁵

Finally, we offer the caveat that our analysis relies on a number of simplifying assumptions. We already mentioned that the synthetic spectra we used were obtained with the Monte Carlo code by Nobili et al. (2008), which was designed to solve radiation transport in a globally twisted magnetosphere. Even though we took thermal photons to originate mostly in a limited polar region, this does not self-consistently describe resonant up-scattering in a magnetosphere where only a limited bundle of field lines is actually twisted, as is probably the case in AXPs (Beloborodov 2009). Moreover, as we discussed in Section 4.2.2, the ntzang XSPEC model is available only in tabular form and it was created assuming emission at constant temperature from the entire star surface. As such, it is not suited to be applied directly to the present case. As a compromise, we decided to fit the spectra by adding together two/three (absorbed) ntzang components, each associated with one of the emitting regions, at temperatures T_h, T_w, and T_c, respectively. While this procedure works (and is routinely employed) in the case of blackbody spectra, it is expected to be only approximately correct when different ntzang components are added together. The reason is that the effects of resonant scattering on thermal photons depend on the location of the primary emission, since the magnetospheric electron density is not isotropic. As a consequence, assuming thermal emission from a cap of limited size or from the entire star, even if the two are taken at the same temperature, will give rise to different spectra. We checked this approximation for all the spectra we analyzed, finding that the maximum relative error is ∼0.6, while the energy-averaged error is always between 0.2 and 0.4 both for XTE J1810−197 and CXOU J164710.2−455216. An example is shown in Figure 12. Although we are aware that this is not optimal, it provides a reasonable way to describe radiation coming from a magnetar with non-uniform thermal emission within the context of our model.

Zoom In Zoom Out Reset image size

The monitoring of the two TAXPs XTE J1810−197 and CXOU J164710.2−455216, carried out with XMM-Newton in recent years, gave us the possibility to test the twisted magnetosphere model and understand how the physical parameters in the two sources change during the post-outburst evolution. We summarize our main findings below, remarking again that they were obtained under a number of assumptions (e.g., globally twisted field, three temperature thermal maps).

• 1.  
  Soon after the outburst onset the surface thermal distribution in XTE J1810−197 and CXOU J164710.2−455216 is well described by three components: a hot cap and a surrounding warm corona while the rest of the neutron star surface is at a lower temperature.
• 2.  
  The analysis of the pulse profile evolution for XTE J1810−197 revealed that both the hot cap and the warm corona decrease in size so that in the last observation (2007 September) virtually all the neutron star surface emits at a temperature compatible with the quiescent one.
• 3.  
  The same analysis for CXOU J164710.2−455216 showed that the hot cap decreases in temperature and size, while the warm corona remains constant in temperature while it increases in size. In the last two observations we examined (2008 August and 2009 August), the source thermal map comprises a hot cap covering ∼30% of the neutron star surface, while the remaining surface is cooler. There are hints that this could be the quiescent state of the TAXP.
• 4.  
  For both sources the twist angle is highest at the outburst onset and then monotonically decreases in time until it reaches a nearly constant, non-zero value.
• 5.  
  The same model configuration which best-fits the observed pulse profiles (thermal map, twist angle, electron bulk velocity, and geometrical angles) provides a reasonable description of XMM-Newton spectra in the 0.1–10 keV band for both sources.

To our knowledge this is the first time that a self-consistent spectral and timing analysis, based on a realistic modeling of resonant scattering, was carried out for magnetar sources, considering simultaneously a large number of data sets over a baseline of years. Present results support a picture in which only a limited portion of the magnetosphere was affected by the twist. Future developments will require detailed spectral calculations in a magnetosphere with a localized twist which decays in time.

We are grateful to an anonymous referee for his/her constructive criticism and helpful suggestions which helped in improving a previous version of this paper. This work is partially supported by INAF-ASI through grant AAE I/088/06/0.