Unveiling the thermal and magnetic map of neutron star surfaces though their X-ray emission: method and light-curve analysis

Recent Chandra and XMM–Newton observations of a number of X-ray ‘dim’ pulsating neutron stars have revealed quite unexpected features in the emission from these sources. Their soft thermal spectrum, believed to originate directly from the star surface, shows evidence for a phase-varying absorption line at some hundred eVs. The pulse modulation is relatively large (pulsed fractions in the range ∼ 12–35 per cent), the pulse shape is often non-sinusoidal, and the hard X-ray colour appears to be anticorrelated in phase with the total emission. Moreover, the prototype of this class, RX J0720.4 − 3125, has been found to undergo rather sensible changes in both its spectral and timing properties over a time-scale of a few years. All these new ﬁndings seem difﬁcult to reconcile with the standard picture of a cooling neutron star endowed with a purely dipolar magnetic ﬁeld, at least if surface emission is produced in an atmosphere on top of the crust. In this paper we explore how a dipolar + quadrupolar star-centred ﬁeld inﬂuences the properties of the observed light curves. The phase-resolved spectrum has been evaluated accounting for both radiative transfer in a magnetized atmosphere and general relativistic ray-bending. We computed over 78 000 light curves, varying the quadrupolar components and the viewing geometry. A comparison of the data with our model indicates that higher-order multipoles are required to reproduce the observations.


I N T RO D U C T I O N
Over the last few years, a number of high-resolution spectral and timing observations of thermally emitting neutron stars (NSs) have become available thanks to new-generation X-ray satellites (both Chandra and XMM-Newton), opening new perspectives in the study of these sources. Thermal emission from isolated NSs is presently observed in more than 20 sources, including active radio pulsars, soft γ -repeaters, anomalous X-ray pulsars, Geminga and Gemingalike objects, and X-ray dim radio-silent NSs. There is by now a wide consensus that the soft, thermal component directly originates from the surface layers as the star cools down. If properly exploited, the information it conveys is bound to reveal much about the physics of NSs, shedding light on their thermal and magnetic surface distribution and ultimately probing the equation of state of matter at supra-nuclear densities.
Although thermal surface emission seems indeed to be a ubiquitous feature in isolated NSs, a power-law, non-thermal compo-E-mail: sz@mssl.ucl.ac.uk (SZ); turolla@pd.infn.it (RT) nent (likely produced in the star magnetosphere) is present in most sources, where it often dominates the X-ray spectrum. Moreover, the intrinsic X-ray emission from young radio pulsars may be significantly contaminated by the contribution of the surrounding supernova remnant. In this respect, the seven X-ray dim isolated neutron stars (XDINSs) discovered by ROSAT are a most notable exception. In a sense, we can claim that these are the only 'genuinely isolated' NSs and their soft thermal emission is unmarred by (non-thermal) magnetospheric activity or by the presence of a supernova remnant or a binary companion (see, for example, Treves et al. 2000 andHaberl 2004 for reviews;Zane et al. 2005). XDINSs play a key role in compact objects astrophysics: these are the only sources in which we can have a clean view of the compact star surface, and as such offer an unprecedented opportunity to confront theoretical models of NS surface emission with observations. The X-ray spectrum of XDINSs is with no exception blackbodylike with temperatures in the range ∼40-100 eV and, thus far, pulsations have been detected in five sources, with periods in the range 3-11 s (see Table 1 and references therein). In each of the five cases, the pulsed fraction is relatively large (∼12-35 per cent). Quite surprisingly, and contrary to what we would expect in a simple dipolar Table 1. Isolated NS parameters. The light curves used in this paper are taken from the references listed in the fifth column, which correspond to: 1, ; 2, Haberl et al. (2003);3, de Vries et al. (2004);4, Zane et al. (2005 geometry, often the hardness ratio is minimum at the pulse maximum (Cropper et al. 2001;Haberl et al. 2003). Broad absorption features have been detected around ∼300-700 eV in all pulsating XDINSs and the line strength appears to vary with the pulse phase.
In addition, the X-ray light curves exhibit a certain asymmetry, with marked deviations from a pure sinusoidal shape at least in the case of RBS 1223 (Haberl et al. 2003;Schwope et al. 2005). XDINSs were unanimously believed to be steady sources, as indicated by several years of observations for the brightest of them. Unexpectedly, and for the first ever time, XMM-Newton observations have recently revealed a substantial change in the spectral shape and pulse profile of the second most luminous source, RX J0720.4−3125, over a time-scale of ∼2 yr Vink et al. 2004). Possible variations in the pulse profile of RX J0420.0−5022 over a similar time-scale (∼0.5 yr) have also been reported, although only at a low significance level ).
In the standard picture, emission from an isolated, cooling NS arises when thermal radiation originating in the outermost surface layers traverses the atmosphere which covers the star crust. Although the emerging spectrum is thermal, it is not a blackbody because of radiative transfer in the magnetized atmosphere and the inhomogeneous surface temperature distribution. The latter is controlled by the crustal magnetic field, because thermal conductivity across the field is highly suppressed, and bears the imprint of the field topology. Besides the spectrum, radiative transfer and the surface temperature distribution act together in shaping the X-ray light curve. Pulse profiles produced by the thermal surface distribution induced by a simple core-centred dipolar magnetic field were investigated long ago by Page (1995), under the assumption that each surface patch emits (isotropic) blackbody radiation. Because of gravitational effects and because of the smooth temperature distribution (the temperature monotonically decreases from the poles to the equator), the pulse modulation is quite modest (pulsed fraction 10 per cent) for reasonable values of the star radius. Moreover, with the temperature distribution being symmetric about the magnetic equator, the pulse shape itself is always symmetrical, regardless of the viewing geometry. Larger pulsed fractions may be reached by the proper inclusion of an atmosphere. In fact, in a strongly magnetized medium photon propagation is anisotropic and occurs preferentially along the field (magnetic beaming; e.g. Pavlov et al. 1994). Nevertheless, retaining a dipolar temperature distribution will always result in a symmetric pulse profile.
The quite large pulsed fraction, pulse asymmetry, and possibly long-term variations, recently observed in XDINSs seem therefore difficult to explain by assuming that the thermal emission originates at the NS surface, at least when assuming that the thermal surface distribution is that induced by a simple core-centred dipolar magnetic field. It should be stressed that, although the dipole field is a convenient approximation, the real structure of the NS magnetic field is far from being understood; for example, it is still unclear if the field threads the entire star or is confined in the crust only (e.g. Geppert, Küker & Page 2004 and references therein). Whatever the case, there are both observational and theoretical indications that the NS surface field is 'patchy' (e.g. Geppert, Rheinhardt & Gil 2003;Urpin & Gil 2004 and references therein). The effects of a more complex field geometry have been investigated by Page & Sarmiento (1996), who considered a star-centred dipole+quadrupole field, again assuming isotropic emission. The presence of multipolar components induces large temperature variations even between nearby regions, and this results in larger pulsed fractions and asymmetric pulse profiles.
The high-quality data now available for thermally emitting NSs, and XDINSs in particular, demand for a detailed modelling of surface emission to be exploited to a full extent. Such a treatment should combine both an accurate formulation of radiation transport in the magnetized atmosphere and a quite general description of the thermal and magnetic surface distributions, which, necessarily, must go beyond the simple dipole approximation. The ultimate goal is to produce a completely self-consistent model, capable of reproducing simultaneously both the spectral and timing properties. In this paper we take a first step in this direction and present a systematic study of X-ray light curves from cooling NSs, accounting for both a quadrupolar magnetic field (in addition to the core-centred dipole) and radiative transfer in the magnetized atmosphere. We computed over 78 000 model light curves, exploring the entire parameter space, in both the geometrical angles and the quadrupolar components. This large data set has been analysed using multivariate statistical methods (principal component analysis and cluster analysis) and we show that a non-vanishing quadrupolar field is required to reproduce the observed XDINS pulse profiles.

Going quadrupolar
In this section we describe the approach we use to compute the phase-dependent spectrum emitted by a cooling NS as seen by a distant observer. This issue has been addressed by several authors in the past under different assumptions, and basically divides into two steps. The first involves the computation of the local (i.e. evaluated by an observer at the star surface) spectrum emitted by each patch of the star surface, while the second requires the collection of the contributions of surface elements which are 'in view' at different rotation phases, taking proper account of the fact that only rays propagating parallel to the line of sight actually reach the distant observer. Details on each are presented in Sections 2.2 and 2.3, and further in Appendix A; here we discuss some general assumptions which are at the basis of our model.
We take the NS to be spherical (mass M and radius R) and rotating with constant angular velocity ω = 2π/P, where P is the period. Because XDINSs are slow rotators (P ≈ 10 s), we can describe the space-time outside the NS in terms of the Schwarzschild geometry (see, for example, Cadeau, Leahy & Morsink 2005 for a more complete discussion about the effects of rotation).
The star magnetic field is assumed to possess a core-centred dipole+quadrupole topology, B =B dip + B quad . Introducing a polar coordinate system whose axis coincides with the dipole axis, the (polar) components of the dipole field at the star surface are where B p is the field strength at the magnetic pole, and θ and φ are the magnetic colatitude and azimuth. The functions f dip and g dip account for the effect of gravity and depend on the dimensionless star radius x ≡ R/R S with R S = 2GM/c 2 ; their explicit expressions can be found in Page & Sarmiento (1996, see also references therein). The quadrupolar field can easily be derived from the spherical harmonics expansion and its expression, again at r = R, can be cast in the form where q i are arbitrary constants. The polar components of the five generating vectors B (i) quad are reported in Page & Sarmiento (1996). We just note that their expression for the radial component of the zeroth vector contains a misprint and should read B (0) quad,r = (3 cos 2 θ − 1)/2. General relativistic effects are included in the quadrupolar field by multiplying the radial and angular components by the two functions f quad (x) and g quad (x) respectively (see again Page & Sarmiento 1996 for their expressions and further details).
The NS surface temperature distribution, T s , will in general depend on how heat is transported through the star envelope. Under the assumption that the field does not change much over the envelope scaleheight, heat transport can be treated (locally) as onedimensional. The surface temperature then depends only on the angle between the field and the radial direction, cos α =B · n, and on the local field strength B (see Page 1995). As shown by Greenstein & Hartke (1983), a useful approximation is to write where the ratio of the conductivities perpendicular (K ⊥ ) and parallel (K ) to the field is assumed to be constant. The polar value T p fixes the absolute scale of the temperature and is a model parameter (Page 1995 and references therein). For field strengths 10 11 G, the conductivity ratio is much less than unity and equation (5) simplifies to This expression is used in the present investigation. An example of the thermal surface distribution induced by a quadrupolar field is shown in Fig. 1. Different approaches which account for the variation of the conductivities (e.g. Heyl & Hernquist 1998) yield similar results. Quite recently, Geppert et al. (2004) have investigated the influence of different magnetic field configurations on the surface temperature distribution. They have found that, contrary to starcentred core fields, crustal fields may produce steeper surface temperature gradients. The inclusion of this effect is beyond the purpose of this first paper. However, we caveat that, for temperatures expected in XDINSs (≈10 6 K), the differences between the two magnetic configurations start to be significant when the field strength is >10 13 G.

Radiative transfer
The properties of the radiation spectrum emitted by highly magnetized, cooling NSs have been thoroughly discussed in the literature (e.g. Shibanov et al. 1992;Pavlov et al. 1994;Zane et al. 2001;Ho & Lai 2001. Because the pressure scaleheight is much smaller than the star radius, model atmospheres are usually computed in the plane parallel approximation. Besides surface gravity, the spectrum emerging from each plane parallel slab depends on both the surface temperature T s and magnetic field B, either its strength or orientation with respect to the local normal, which coincides with the unit vector in the radial direction n. In order to proceed, we introduce a (θ, φ) mesh which naturally divides the star surface into a given number of patches. Once the magnetic field has been specified, each surface element is characterized by a precise value of B, α and T s . The atmospheric structure and radiative transfer can then be computed locally by approximating each atmospheric patch with a plane parallel slab, infinitely extended in the transverse direction and emitting a total flux σ T 4 s .
Radiative transfer is solved numerically using the approach described in . The code is based on the normal mode approximation for the radiation field propagating in a strongly magnetized medium and incorporates all relevant radiative processes. The full angle and energy dependence in both the plasma opacities and the radiation intensity is retained. In this respect we note that for an oblique field [i.e. (B/B) ·n = 1] the intensity is not symmetric around n and depends explicitly on both propagation angles. If k is a unit vector along the photon direction, at depth τ in the atmosphere the intensity has the form I = I E (τ , μ, ϕ) where E is the photon energy, μ = n · k and ϕ is the associated azimuth. Calculations are restricted to a completely ionized H atmosphere (see  for a treatment of ionization in H atmospheres).
Because the numerical evaluation of model atmospheres is computationally quite demanding, especially for relatively high, oblique fields and T s < 10 6 K, we preferred to create an archive beforehand, by computing models for pre-assigned values of cos α, B and T s . The range of the latter two parameters should be wide enough to cover the surface variation of B and T in all the cases of interest: 12 log B 13.5 and 5.4 log T s 6.6. Particular care was taken to generate models in the entire α domain, 0 cos α 1. According to the adopted surface temperature distribution (equation 6), the regions close to the equator have very low temperatures and cannot be associated with any model in the archive. However, being so cool, their contribution to the observed spectrum is negligible (see Section 2.4). For each model, the emerging radiation intensity, I E (τ = 0, μ, ϕ), is stored, for 0.01 E 10 keV, 0 μ 1 and 0 ϕ 2π. The archive itself consists of a sixdimensional array I(B, T s , cos α; E, μ, ϕ) which associates at each set of the parameters B, cos α, T s the (discrete) values of the angleand energy-dependent intensity. Actually, because the code makes use of adaptive computational meshes, the emerging intensities have been interpolated on common energy and angular grids before storage. The final array contains the emerging intensity for about 100 atmospheric models, evaluated at 100 energy bins and on a 20 × 20 angular mesh, (μ, ϕ).

Observed spectrum
The problem of computing the pulse profile produced by hot caps on to the surface of a slowly rotating NS including gravity effects was first tackled by Pechenick, Ftaclas & Cohen (1983). Their results were then generalized to the case of emission from the entire star surface with an assigned temperature distribution by Page (1995) and Lloyd et al. (2004), in the case of isotropic and non-isotropic radiation fields, respectively. The approach used here follows that discussed in the two papers quoted above, which we refer to for more details. For the sake of clarity, we present a Newtonian derivation first. Relativistic ray-bending can then be accounted for quite straightforwardly.
The NS viewing geometry is described in terms of two angles χ and ξ , which give the inclination of the line of sight and of the dipole axis with respect to the star spin axis. Let z, b dip and p denote the unit vectors along the same three directions. Let us moreover introduce two Cartesian coordinate systems, both with origin at the star centre: the first, (X, Y, Z), is fixed and such that the Z-axis coincides with the line of sight while the X-axis is in the (z, p) plane; the second, (x, y, z), rotates with the star. The z-axis is parallel to b dip while the choice of the x-axis will be made shortly. Each Cartesian frame has an associated polar coordinate system, with polar axes along Z and z, respectively. The colatitude and azimuth are ( , ) in the fixed frame, and ( , ) in the rotating one (the latter are just the magnetic colatitude and azimuth introduced in Section 2.1).
In the following we express vectors through their components: these are always the Cartesian components referred to the fixed frame, unless otherwise explicitly stated. The same components are used to evaluate both scalar and vector products. Upon introducing the phase angle γ = ωt, it follows from elementary geometrical considerations that p = (−sin χ , 0, cos χ) and b dip = (−sin χ cos ξ − cos χ sin ξ cos γ , − sin ξ sin γ , cos χ cos ξ + sin χ sin ξ cos γ ). It can easily be verified that q = (cos ξ sin γ , cos γ , sin ξ cos γ ) represents a unit vector orthogonal to p and rotating with angular velocity ω. We then choose the x-axis in the direction of the (normalized) vector component of q perpendicular to b dip the y-axis is parallel to b dip × q ⊥ . The local unit normal relative to a point on the star surface of coordinates ( , ) is readily expressed as n = (sin cos , sin sin , cos ). By introducing the unit vector n ⊥ , defined in strict analogy with q ⊥ (see equation 7), the two expressions provide the relations between the two pairs of polar angles, the geometrical angles ξ , χ and the phase. While direct inversion of equation (8) poses no problems because it is 0 θ π, care must be taken to ensure that φ, as obtained from equation (9), covers the entire range [0, 2π]. This is achieved by replacing We are now in the position to compute the total monochromatic flux emitted by the star and received by a distant observer. This is done by integrating the specific intensity over the visible part of the star surface at any given phase (see, for example, Lloyd et al. 2004) where u = sin . Further integration of equation (10) over γ provides the phase-averaged spectrum. As discussed in Section 2.2, the intensity depends on the properties of the surface patch and on the photon direction. The magnetic field strength B can be directly computed from the polar components of B (see Section 2.1 and Page & Sarmiento 1996). The magnetic tilt angle α and the surface temperature (see equation 6) follow from cos α = n · B/B = B r /B, where n is the unit radial vector. The local values of B and T s depend on (θ, φ). They can then easily be expressed in terms of ( , ) for any given phase using equations (8) and (9). Because the star appears point-like, there is a unique ray which reaches the observer from any given surface element and this implies that also μ and ϕ are functions of ( , ). Clearly μ = cos , while ϕ is given by cos ϕ =m · v. The two unit vectors which enter the previous expression are, respectively, the projections of B and z on the plane locally tangent to the surface. They are expressed as m = (cos cos , cos sin , −sin ) and v = (B/B − n cos α)/sin α. The Cartesian components of B needed to evaluate cos ϕ are derived in Appendix A. Gravity effects (i.e. relativistic ray-bending) can now be included in a very simple way. The local value of the colatitude is, in fact, related to that measured by an observer at infinity by the 'ray-tracing' where x = R/R S . Because we are collecting the contributions of all surface elements seen by a distant observer, each patch is labelled by the two angles¯ and . This means that the integrand in equation (10) is to be evaluated precisely at the same two angles and is tantamount to replacing with¯ in all previous expressions. Note, however, that the innermost integral in equation (10) is always computed over . Effects of radiative beaming are illustrated in Fig. 2, where we compare phase-averaged spectra and light curves computed by using radiative atmosphere model with those obtained under the assumption of isotropic blackbody emission. As we can see, the pulse profiles are substantially different in the two cases. Also, accounting for radiative beaming allows us to reach relatively large pulse fractions (∼20 per cent).

Numerical implementation
The numerical evaluation of the phase-dependent spectrum has been carried out using an IDL script. Because our immediate goal is to check if the observed pulse profiles can be reproduced by our model, we start by computing the light curve in a given energy band, accounting for interstellar absorption and the detector response function. Most of the recent observations of X-ray emitting INSs have been obtained by XMM-Newton, so here we refer to the EPIC-pn response function. Both absorption and the detector response depend upon the arrival photon energyĒ = E √ 1 − 1/x, so the pulse profile in the [Ē 1 ,Ē 2 ] range is given by where A is the response function, N H is the column density, and σ is the interstellar absorption cross-section (e.g. Morrison & McCammon 1983). Because the energy integral J does not involve geometry, it is evaluated first. We select three energy bands, corresponding to the soft (0.1-0.5 keV) and hard (0.5-1 keV) X-ray colours, and to the total range 0.1-5 keV. Results are stored as described in Section 2.2 in the case of the quantity I, the only difference being that energy index in the array J runs now from 1 to 3 corresponding to the three energy intervals defined above. We then introduce a (μ = cos , ) mesh by specifying 50 × 50 equally spaced points in the [0, 1] and [0, 2π] intervals, and we interpolate the intensity array at the required values of μ, μ = u. Next, the values of¯ relative to the u grid are computed from equation (11). All these steps can be performed in advance and once for all, because they do not depend on the viewing geometry or the magnetic field.
Then, once the angles χ , ξ and the phase γ have been specified, the magnetic colatitude and azimuth, θ (¯ , , γ ), φ(¯ , , γ ) can be evaluated. We use 32 phase bins and a set of five values for each of the two angles χ , ξ = (0 • , 30 • , 50 • , 70 • , 90 • ). The magnetic field is assigned by prescribing the strength of the quadrupolar components relative to the polar dipole field, b i = q i /B p (i = 1, . . . , 5), in addition to B p itself; in our grid, each b i can take the values (0, ± 0.25, ± 0.5). For each pair θ , φ we then compute B and cos ϕ; cos α gives the surface temperature T s once T p has been chosen.
The corresponding values of the intensity are obtained by linear interpolation of the array J . Surface elements emitting a flux two orders of magnitudes lower than that of the polar region (σ T 4 p ) were assumed not to give any contribution to the observed spectrum. Finally, direct numerical evaluation of the two angular integrals in equation (12) gives the light curves. Although in view of future application we computed and stored the light curves in the three energy bands mentioned above, the results presented in Sections 3 and 4 are obtained using always the total (0.1-5 keV) energy band. To summarize, each model light curve depends on χ , ξ , and the five b i . No attempt has been made here to vary also T p and B p , which have been fixed at 140 eV and 6 × 10 12 G, respectively. A total of 78 125 models have been computed and stored. Their analysis is discussed in Section 3.

A NA LY S I N G L I G H T C U RV E S A S A P O P U L AT I O N
As discussed in Section 2, under our assumptions the computed light curve is a multidimensional function which depends in a complex way on several parameters. Therefore, an obvious question is whether or not we can identify some possible combinations of the independent parameters that are associated with particular features observed in the pulse shape. The problem of quantifying the degree of variance of a sample of individuals (in our case the light curves) and identifying groups of 'similar' individuals within a population is a common one in behavioural and social sciences. Several techniques have been extensively detailed in many books of multivariate statistics (e.g. Kendall 1957;Manly 1998) and, although they have been little used in physical sciences, a few applications to astrophysical problems have been presented over the past decades (see, for example, Whitney 1983; Mittaz, Penston & Snijders 1990;Heyer & Schloerb 1997).
We focus here on a particular tool called principal component analysis (PCA), which appears promising for a quantitative classification of the light-curve features. The main goal of PCA is to reduce the number of variables that need to be considered in order to describe the data set, by introducing a new set of variables z p -called the principal components (PCs) -which can discriminate most effectively among the individuals in the sample (i.e. the light curves in our present case). The PCs are uncorrelated and mutually orthogonal linear combinations of the original variables. Besides, the indices of the PCs are ordered so that z 1 displays the largest amount of variation, z 2 displays the second largest, and so on; that is, var (z 1 ) var (z 2 ) . . . var (z p ) where var (z k ) is the variance in the sample associated with the kth PC. Although the physical meaning of the new variables may not, in general, be immediate, it often turns out that a good representation of the population is obtained by using a limited number of PCs, which allows us to treat the data in terms of their true dimensionality. Beyond that, PCA represents a first step towards other types of multivariate analyses, such as cluster analysis. This is a tool which allows the identification of subgroups of objects so that 'similar' ones belong to the same group. When applying the cluster analysis algorithm, PCA is performed first in order to reduce the number of original variables down to a smaller number of PCs. This can substantially reduce the computational time.
Because both tools are extensively described in the literature, we omit all mathematical details, for which an interested reader may refer, for example, to Kendall (1957) and Manly (1998). Let us denote by y ps ( p = 1, . . . , P; s = 1, . . . , S) the values of the intensity computed at each phase bin, for each model light curve. Let us also introduce the 'centred' data x ps , defined as where μ p and s p are the mean and the variance of the computed data, respectively. In order to shift to the new PC variables z ps , we computed the transformation matrix V such that A sufficient condition to specify univocally V is to impose that the axes of the new coordinate system (i.e. the PCs) are mutually orthogonal and linearly independent. By applying PCA to our sample of models, we have found that each light curve can be reproduced using only the first 20-21 more significant PCs (instead of 32 phases) and that the first four PCs alone account for as much as ∼85 per cent of the total variance. Moreover, ∼72 per cent of the variance is in the first three PCs only. It is therefore meaningful to introduce a graphical representation of the model data set in terms of the first three z i . This is shown in Fig. 3 where black/red squares give the position in the z 1 z 2 z 3 space of quadrupolar/dipolar models. To better visualize the latter, an additional set of light curves was computed, bringing the total number of dipolar models displayed in Fig. 3 to 100.
An insight into the light-curve property of each of the PC measures can be obtained by inspecting the coefficients v pq of the linear combination, which gives z p for an assigned data set (see equation 14). Because z p = q v pq x q ∝ 2π 0 v p (γ )F(γ ) dγ , this is tantamount to assessing the effect of the filter v p (γ ) on the light curve F(γ ) knowing the values of the former at a discrete set of phases. The first four v p are shown in Fig. 4. The first PC provides a measure of the amplitude of the pulse; it is always z 1 > 0, and large values of z 1 correspond to low pulsed fractions. Both z 2 and z 3 may take either sign (same for higher-order PCs) and give information on the pulse shape. We note that, although the absolute phase is in principle an arbitrary quantity, the whole model population has been computed using the same value. Therefore, when studying the morphological properties of the sample of light curves, it is meaningful to refer to the symmetry properties with respect to the parameter γ . Large and negative values of z 2 imply that the main contributions to the light curve come from phases around zero. z 3 measures the parity of the light curve with respect to the half period. Pulses which are symmetric have z 3 = 0. As Fig. 3 clearly shows, the PCA is very effective in discriminating purely dipolar from quadrupolar models. The former cluster in the 'tip of the cloud', i.e. at large values of z 1 , negative values of z 2 and in the z 3 = 0 plane, as expected. In fact, dipolar pulse patterns are always symmetrical and their pulsed fraction is quite small (semi-amplitude 10 per cent). It is worth noticing that quadrupolar magnetic configurations too can produce quite symmetrical light curves, e.g. the black squares in Fig. 3 with z 3 ∼ 0. However, in this case the pulsed fraction may be much  larger, as indicated by the lower values of z 1 attained by quadrupolar models with z 3 = 0 in comparison with purely dipolar models. This implies that a symmetrical light curve is not per se an indicator of a dipolar magnetic geometry.
As suggested by some authors (e.g. Heck 1976;Whitney 1983), PCs can then be used as new variables to describe the original data. However, in the case at hand, the problem is that although PCs effectively distinguish among pulse patterns, they have a purely geometrical meaning and cannot be directly related to the physical parameters of the model (b i , χ, ξ ). We have found that the standard regression method does not allow us to link the PCs with the model parameters, which is likely to be a signature of a strong nonlinearity (see Section 5). Instead, the PCA can be regarded as a method to provide a link between the numerous light curves, in the sense that models 'close' to each other in the PC space will have similar characteristics. Unfortunately, despite the fact that different definitions of metrics have been attempted, so far we have found it difficult to translate the concept of 'proximity' in the PC space into a corresponding 'proximity' in the seven-dimensional space of the physical parameters ξ , χ and b i (i = 1, . . ., 5). By performing cluster analysis, we have found that two separate subgroups are clearly evident in the PC space, one of which encompasses the region occupied by purely dipolar models (see Fig. 5, top panel). However, again it is not immediate to find a corresponding subdivision in the physical space. Because of these difficulties, we postpone a more detailed statistical analysis to a follow-up paper, and, as discussed in the next section, we concentrate on the direct application of our models to the observed light curves of some isolated NSs.

A N A P P L I C AT I O N TO X -R AY D I M I S O L AT E D N E U T RO N S TA R S
In the light of the discussion at the end of the previous section, the PCA may be used to gain insights into the ability of the present model to reproduce observed light curves. A simple check consists of deriving the PC representation of the pulse profiles of a number of sources and verifying if the corresponding points in the z 1 z 2 z 3 space fall inside the volume covered by the models. We stress once again that, although the model light curves depend upon all the PCs, z 1 , z 2 and z 3 alone provide a quite accurate description of the data set because they account for a large fraction of the variance. In this sense the plots in Fig. 3 give a faithful representation of the light curves, with no substantial loss of information: profiles close to each other in this three-dimensional space exhibit a similar shape. To this end, we took the published data for the first four pulsating XDINSs listed in Table 1 and rebinned the light curves to the same 32 phase intervals used for the models. Because the PCA of the model data set already provided the matrix of coefficients v pq (see equation 14), the PCs of any given observed light curve are z obs p = q v pq x obs q , where x obs q is the (standardized) X-ray intensity at phase γ q . As can be seen from Fig. 3, the observed pulse profiles never fall in the much smaller region occupied by purely dipolar models. However, they all lie close to the quadrupolar models, indicating that a quadrupolar configuration able to reproduce the observed features may exist. A possible exception to this is the first observation of RX J0720.4−3125 (XMM-Newton revolution 078), for which the pulse profile appears quite symmetric and the pulsed fraction is relatively small (see Table 1). While a purely dipolar configuration may be able to reproduce the light curve for revolution 078, a visual inspection of Fig. 3 shows that this is not the case for the second observation of the same source (XMM-Newton revolution 711). Despite the fact that the pulse shape is still quite symmetrical, the pulsed fraction is definitely too large to be explained, within the current model, with a simple dipole field. The same considerations apply to the light curve of RX J0420.0−5022.
Just as a counterexample, we added in Fig. 3 the PC representation of the X-ray light curve of the anomalous X-ray pulsar 1E 1048.1−5937 observed in two different epochs, 2000 December and 2003 June (see Mereghetti et al. 2004). The new data points fall well outside the region populated by quadrupolar models and can be seen only in the last panel of Fig. 3. For this source, the pulsed fraction is so high (89 and 53 per cent in June and December, respectively) that no quadrupolar configuration could account for it. In terms of PCs, both points have a relatively low value of z 1 (z 1 = 9.9 and z 1 = 6.9). In this case no fit can be found, not surprisingly, because we do expect a large contribution from a non-thermal power-law component to the X-ray emission of anomalous X-ray pulsars.
To better understand to what extent quadrupolar surface distributions may indeed produce light curves close to those observed, we select the model in the data set which is 'closer' to each of the observed light curves. This is done by looking first for the minimum Euclidean distance between the observed and the model pulse profiles in the PC space. Note that in this case all the relevant PCs were used; that is to say, the distance is defined by D 2 = 20 p=1 (z p −z obs p ) 2 . The computed model which minimizes D is then used as the starting point for a fit to the observed light curve which yields the best estimate of the model parameters. The fitting is performed by computing 'on the fly' the light curve for the current values of the parameters and using the standard (i.e. not the PC) representation of both model and data. The quadrupolar components and viewing angles The results of the fits are shown in Figs 5, 6 and 7 for B dip = 6 × 10 12 G, log T p (K) = 6.1-6.2. It is apparent that the broad characteristics of all the XDINS light curves observed so far may be successfully reproduced for a suitable combination of quadrupolar magnetic field components and viewing angles. However, although in all cases a fit exists, we find that in general it is not necessarily unique. This means that the model has no 'predictive' power in providing the exact values of the magnetic field components and viewing angles. For this reason we do not attempt a more complete fit, i.e. leaving also T p and B dip free to vary, nor do we derive parameter uncertainties or confidence levels. Our goal at this stage has been to show that there exists at least one (and more probably more) combination(s) of the parameters that can explain the observed pulse shapes, while this is not possible assuming a pure dipole configuration.
The case of RX J0720.4−3125 deserves, however, some further discussion. This source, which was believed to be stationary and as such included among XMM-Newton EPIC and RGS calibration sources, was recently found to exhibit rather sensible variations in both its spectral and timing properties Vink et al. 2004). In particular, the pulse shape changed over the ∼2 yr time between XMM revolutions 78 and 711. de  proposed that the evolution of RX J0720.4−3125 may be produced by a (freely) precessing NS. This scenario can be tested by our analysis, because in a precessing NS only the two angles ξ and χ are expected to vary while the magnetic field remains fixed. This means that having found one combination of parameters which fits the light curve of revolution 78, a satisfactory fit for revolution 711 should be obtained for the same b i and different values of the angles. Despite several attempts, in which the proper values of T p as derived from the spectral analysis of the two observations were used, we were unable to obtain a good fit varying the angles only (see Fig. 8). We also performed a simultaneous fit to both light curves starting from a general trial solution and asking that b i are the same (but not necessarily coincident with those that best fit data from revolution 78), while the two angles (and the initial phases)  Fig. 7, lower panel. In this particular example, we also allowed a fractional change of ∼20 per cent in the temperature between the two revolutions. Both attempts give an unsatisfactory result.
are kept distinct (see Fig. 9). Both approaches clearly indicate that a change in both sets of quantities (magnetic field and angles) is required (as in Fig. 7). From a physical point of view it is not clear how magnetic field variations on such a short time-scale may be produced; therefore, at present no definite conclusion can be drawn. For instance, one possibility that makes conceivable a change of the magnetic field structure and strength on a time-scale of years is that the surface field is small scaled (Geppert et al. 2003). In this case, even small changes in the inclination between the line of sight and local magnetic field axis may cause significant differences in the 'observed' field strength.

D I S C U S S I O N
XDINSs may indeed represent the Rosetta stone for understanding many physical properties of NSs at large, including their equation of state. Despite their potential importance, only recently have detailed observations of these sources become progressively available with the advent of the Chandra and XMM-Newton satellites. These new data, while confirming the thermal, blackbody-like emission from the cooling star surface, have revealed a number of spectral and timing features, which have opened a new window on the study of these objects. Some issues are of particular importance in this respect: (i) the discovery of broad absorption features at a few hundred eVs; (ii) the quite large pulsed fractions; (iii) the departure of the pulse shape from a sine wave; (iv) the long-term evolution of both the spectral and timing properties seen in RX J0720.4−3125 and, to some extent, in RX J0420.0−5022. Pulse-phase spectroscopy confirms that spectral and timing properties are interwoven in a way which appears more complex than that expected from a simple misaligned rotating dipole, as the anticorrelation of the absorption-line strength and of the hardness ratio with the intensity along the pulse testify.
Motivated by this, we have undertaken a project aimed at studying systematically the coupled effects of (i) radiative transfer in a magnetized atmospheric layer, (ii) thermal surface gradients, and (iii) different topologies of the surface magnetic field in shaping the spectrum and pulse shape. The ultimate goal of our investigation is to obtain a simultaneous fit of both pulse profile and (phase-averaged and phase-resolved) spectral distribution. As detailed comparisons of synthetic spectra with observations have shown, no completely satisfactory treatment of spectral modelling for these sources is available as yet. For this reason, here we have presented the general method and concentrated on the study of the light curves, computed assuming a pure H, magnetized atmosphere. The pulse shapes, in fact, should be less sensitive to the details of the chosen atmospheric model.
We provide a caveat for the reader that our results have been computed under a number of simplifying assumptions. For instance, there are still considerable uncertainties in current modelling of the NS surface thermal emission. First, we just mention here that most of the observed NS spectra cannot be reproduced by the theoretical models currently available (see  for reviews and references therein). Secondly, the surface temperature has been derived using equation (5) which is based on the assumption that the temperature profile is only dictated by the heat transferred from the star interior. While this is expected to be the main mechanism, other effects may significantly contribute. For instance, heating of the star surface may occur because of magnetic field decay, and the polar caps may be reheated by back-flowing particles or by internal friction. Thirdly, we have computed the emissivity by assuming that each atmospheric, plane parallel, patch emits a total flux σ T 4 s . In other words, while the spectral distribution is computed using a proper radiative transfer calculation, we have introduced the further assumption that each slab emits the same total flux as a blackbody radiator. A more consistent approach would be to compute the spectrum emitted by each patch by taking the value of T s and the efficiency of the crust as a radiator (Turolla, Zane & Drake 2004) as the boundary condition at the basis of the atmosphere. Our working assumption avoids the burden of creating a complete grid of models in this parameter, with the disadvantage that the spectral properties of each patch may be slightly approximated.
As far as the application to XDINSs is concerned, the greatest uncertainties arise because in this paper we are not fitting simultaneously the observed spectra and pulse shape. The quadrupolar components and viewing angles are treated as free parameters while the polar values T p and B dip are fixed and, of course, they must fall in the domain spanned by our atmosphere model archive. As stated earlier, albeit we have checked that a fit is possible for different combinations of B dip and T p in the allowed range, in general the archive will not contain the exact values of T inferred from spectral observations of the coldest XDINSs. As long as we make the same assumption that the local spectrum emitted by each patch of the star surface is computed using fully ionized, magnetized hydrogen atmosphere models, we do still expect to reach a good fit by using a temperature of a few tens of eV smaller than those used, albeit with different values of the quadrupolar components. However, partial ionization effects are not included in our work, and bound atoms and molecules can affect the results, changing the radiation properties at relatively low T and high B Potekhin et al. 2004).
Even more crucial is the issue of fitting with a realistic value of the magnetic field. The recent detection of absorption features at ≈300-700 eV in the spectrum of XDINSs and their interpretation in terms of proton cyclotron resonances or bound-bound transitions in H or H-like He, may indicate that these sources possess strong magnetic fields, up to B ∼ 9 × 10 13 G van Kerkwijk et al. 2004;Zane et al. 2005). A (slightly) more direct measurement, based on the spin-down measure, has been obtained so far only in one source (i.e. RX J0720.4−3125; e.g. Cropper et al. 2004;Kaplan & van Kerkwijk 2005). In this case the measurement points at a more moderate field, B ∼ (2-3) × 10 13 G, which is only a few times larger than that used here. However, the possibility that (some) of these objects are ultramagnetized NSs is real. Were this to be confirmed, our model archive should be extended to include higher field values. However, this poses a serious difficulty, because the numerical convergence of model atmospheres is particularly problematic at such large values of B for T 10 6 K. Moreover, as mentioned in Section 2.1, if such high field strengths are associated to crustal (and not to star-centred) fields, the surface temperature gradient is expected to be substantially different from that used in the present investigation (Geppert et al. 2004).
The application presented here makes only use of the properties of the pulse profile in the total energy band, and does not exploit colour and spectral information available for these sources. This worsens the problem of finding a unique representation, a problem which is to some extent intrinsic, because of the multidimensional dependence of the fitting function on the various physical parameters. The aim of our future work is to reduce the degeneracy by refining the bestfitting solutions by using information from the light curves observed in different colour bands and/or from the spectral line variations with spin pulse. Also, a more detailed statistical analysis on the model population based on algorithms more sophisticated than the PCA (Gifi 1990;Saegusa 2005) may shed light on the possible correlation between the physical parameters in the presence of strong nonlinearity and possibly on the meaning of the two subclasses identified through a cluster analysis.