{'fraction_non_alphanumeric': 0.039881241257243995, 'fraction_numerical': 0.037369037083559334, 'mean_word_length': 4.33424699253845, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Timing analysis of the isolated neutron star RX J0720.4-3125

We present a combined analysis of XMM-Newton, Chandra and Rosat observations of the isolated neutron star RXJ0720.4-3125, spanning a total period of \sim 7 years. We develop a maximum likelihood periodogramme for our analysis based on the \Delta C-statistic and the maximum likelihood method, which are appropriate for the treatment of sparse event lists. Our results have been checked"a posteriori"by folding a further BeppoSAX dataset with the period predicted at the time of that observation: the phase is found to be consistent. The study of the spin history and the measure of the spin-down rate is of extreme importance in discriminating between the possible mechanisms suggested for the nature of the X-ray emission. The value of \dot P, here measured for the first time, is \approx 10^{-14} s/s. This value can not be explained in terms of torque from a fossil disk. When interpreted in terms of dipolar losses, it gives a magnetic field of B \approx 10^{13} G, making also implausible that the source is accreting from the underdense surroundings. On the other hand, we also find unlikely that the field decayed from a much larger value (B\approx 10^{15} G, as expected for a magnetar powered by dissipation of a superstrong field) since this scenario predicts a source age of \approx 10^4 yrs, too young to match the observed X-ray luminosity. The observed properties are more compatible with a scenario in which the source is \approx 10^6 yrs old, and its magnetic field has not changed substantially over the lifetime.


INTRODUCTION
RXJ0720.4-3125 is a nearby, isolated neutron star (NS) that was originally discovered with ROSAT during a systematic survey of the galactic plane (Haberl et al., 1997). The source exhibits all the common characteristics of the other six ROSAT NS candidates (hereafter dim NSs, see e.g. Treves et al., 2000 for a review): a blackbody-like, soft spectrum with kT ∼ 80 eV; an exceedingly large X-ray to optical flux ratio; a low X-ray luminosity, LX ≈ 10 30 − 10 31 erg s −1 ; a low column density and no evidence for a binary compan- ion. An optical counterpart has been identified (Motch & Haberl, 1998;Kulkarni & van Kerkwijk, 1998).
The true nature of the mechanism powering the X-ray emission from dim NSs is still unclear. Until a few years ago, these sources were thought to constitute a stand alone class and two mechanisms were proposed for their emission: either accretion from the interstellar medium onto an old neutron star or release of thermal radiation from a younger, cooling object (see e.g. Pavlov et al., 1996;Treves et al., 2000;Motch, 2001 and references therein). More recently, it has been noted that the period of RXJ0720.4-3125 is somewhat unusual and falls in the very narrow range in which the anomalous X-ray pulsars (AXPs) periods cluster (P ∼ 6−12 s, see e.g. Mereghetti, 2000). This hints at a possible evolutionary link between dim NSs, AXPs, and soft gammaray repeaters (SGRs, see e.g. Heyl & Kulkarni, 1998;Heyl & Hernquist, 1998;Colpi et al., 2000;Alpar, 1999Alpar, , 2000Alpar, , 2001. Two kind of "unified" scenarios have been then proposed. In the first one, the common mechanism powering the three classes of objects is dissipation of a decaying, superstrong magnetic field (B ∼ > 10 14 − 10 15 G). In this case dim NSs represent the descendants of SGRs and AXPs, and RXJ0720.4-3125 may be one of the closest old magnetars. Alternatively, all three classes may contain neutron stars with lower (canonical) magnetic field (B ≈ 10 12 G) endowed by a fossil disk (Alpar, Ankay & Yazgan, 2001;Alpar, 2001). In this interpretation, dim NSs would be in the propellor phase and would be the progenitors of AXPs and SGRs, the latter having entered an accretion phase.
In any model involving magnetars, proton cyclotron features are expected to lie in the X-ray range, while electron cyclotron lines are expected in the same band for canonical magnetic field. According to Zane et al. (2001), for surface magnetic fields strengths of ∼ 10 14 − 10 15 G spectra exhibit a distinctive absorption feature at the proton cyclotron energy ∼ 0.63(B/10 14 G) keV. The required resolving power is ≈ 100, therefore the detection of this feature is well within the capabilities of XMM-Newton grating spectrometers. Recently, Paerels et al. (2001) have presented spectra of RXJ0720.4-3125 using XMM-Newton, and found that there is no evidence for absorption lines and edges in the X-ray spectrum. Unless different atmospheric effects may concur in suppressing the feature, when taken straightforward the absence of electron or proton cyclotron resonances in the RGS band excludes a range of average magnetic field strengths, B ≈ (0.3 − 2) × 10 11 G and (0.5 − 2) × 10 14 G.
Based on the same XMM-Newton observation, Cropper et al. (2001) presented the pulse-shape analysis of RXJ0720.4-3125, modelling spin pulse intensity and hardness ratio profiles. By assuming that the source of the flux variation is the changing visibility of the heated magnetic polar caps, they derived an upper limit on the cap size and showed that a polar cap larger than ∼ 60 • − 65 • can be rejected at a confidence level of 90%. Whatever the mechanism, the X-ray emitting region is therefore confined to a relatively small fraction of the star surface. It is worth noticing that such a small hot region can not be explained only in terms of non-uniformity in the distribution of the surface cooling temperature induced by an high magnetic field (Greenstein & Hartke 1983;Page 1995). The temperature gradient which is induced by the B-dependence of the thermal conducivity of the neutron star crust is relatively smooth and, even in the (unrealistic) limit case B → ∞ the associated blackbody luminosity drops only by a factor 1/2 at ∼ 60 • and by an order of magnitude at ∼ 77 • . By performing pulse phase spectroscopy, Cropper et al. (2001) also found that the observed hardness ratio is softest slightly before flux maximum. The same characteristic was later discovered by Perna et al. (2001) in the spectra of AXPs. Cropper et al. (2001) suggested two possible explanations for this effect: their best-fitting model is based on radiation beam-ing, while an alternative one assumes a spatially variable absorbing matter, co-rotating in the magnetosphere. The latter may be indeed the case if the star is propelling matter outward (Alpar, Ankay & Yazgan 2001).
Further information about the nature of this puzzling source can be obtained by the spin history. In magnetars magneto-dipolar radiation will cause rapid spin-down at a rateṖ ≈ 10 −11 (B/10 14 G) 2 /P ss −1 , and it has been the positive detection of a secular spin-down of this order that has suggested the association of AXPs and SRGs with ultramagnetized NSs (Kouveliotou et al., 1998(Kouveliotou et al., , 1999Thompson et al., 2001). The preliminary measure ofṖ published by Haberl et al. (1997) for RXJ0720.4-3125 is uncertain to a considerably large value, and does not allow spin-up and spin-down to be discriminated. Very recently, Hambaryan et al. (2002) presented the first evidence of large spin-down rate in a dim NS, RBS1223. Their measure is partially based on ROSAT data where the detection of the period was not highly significant, and also the value (Ṗ = 1.35 +0.69 −0.67 × 10 −11 s/s) is still compatible with being due to torque from a fossil disk. However, it is certainly worth noticing that, when taken face value and interpreted in terms of dipolar rotational losses, this implies the presence of a superstrong magnetic field, B ≈ (3.5 − 6.5) × 10 14 G.
A similar measure ofṖ for the closest pulsating candidate, RXJ0720.4-3125, is therefore of extreme importance, as well as the accurate tracking of its spin history. This is the first opportunity to investigate the pulse timing in the XMM-Newton data, since at the time of the previous papers (Paerels et al. 2001 the timing correlation files were not sufficiently well determined. Here we present a combined analysis of XMM-Newton, Chandra and Rosat data, spanning a period of ∼ 7 years.

OBSERVATIONS
The log of observations used for this study is given in Table 1.
RXJ0720.4-3125 has been observed twice by XMM-Newton: first on 2000 May 13, during the Calibration/Performance Verification phase (Paerels et al. 2001) and later on 21 November 2001, again for calibrations. For the analysis presented in this paper we have used data from all three EPIC cameras on both epochs. The data were reduced using the XMM-Newton Scientific Analysis System (SAS) version 5.1. Source counts in the range 0.12 to 1.2 keV (PN) and 0.1 to 2.0 keV (MOS) were extracted in an aperture of 40 arcsec for the PN X00a, 30 arcsec for PN X00b, and 30 arcsec for the MOS cameras. Events with patterns 0-4 (single and double events) were selected for MOS, while patterns 0-12 (all the valid PN patterns) were retained for PN. Data arriving during episodes of high background which were experienced in the 2000 May 13 observation were excluded. Also in this observation, some PN event times were noted not to be integer multiple of the frame times, but delayed by 1 second: such events (in total about 34 ks from the 50 ks) were corrected in the event list. Finally, times were converted to Barycentric Dynamical Time (TDB) at solar system barycentre.
The Chandra observations were made over three days starting on 2000 Feb 1 using the HRC-S with LETG. The events were extracted from zero-order images of the source (namely, from circles of 2 arcsec radius centered at the star's position). The dispersed data (1 st and higher orders) are strongly contaminated by background (about 50-60%) and are not useful for this task. On the other hand, the background contamination in the zero-order data is negligible, about 2-3%. Also, the zero-order data are not piled-up: this effect can appear only at larger count rates (∼ 10 5 ct/s).
The Rosat observations were taken over a period of five years in 1993, 1996 and 1998. An earlier report of the analysis of the 1993 and 1996 data is in Haberl et al. (1997). The data were re-extracted and re-analysed using the EXSAS software system to make use of the most recent knowledge of spacecraft clock corrections. The 1998 data were taken shortly before a spacecraft emergency safe-hold, and do not benefit from a clock calibration: timing corrections are therefore extrapolated from the last available calibration. Timings may therefore be incorrect up to several seconds (W. Becker, private communication): this means that the reference of this observation to the rest of the dataset is problematic. However, no significant drifts are expected over the duration of the observation itself, so that the data are still useful in a standalone analysis.
Although the countrate is too low for our main analysis, we have also extracted the SAX LECS observation of RX J0720.4-3125 from the SAX data archive (the source was not visible in the other SAX instruments) for use as an a posteriori check on our derived timing parameters. Here we used the event list from the pipeline, correcting the times to the solar system barycentre using the FTOOL earth2sun. This neglects the light travel time from Earth centre to the satellite in low Earth orbit, ∼ 20 ms, which is negligible for the purposes of our a postiori check.
The major datasets in Table 1 are from the two XMM Newton observations which have high count-rates (particularly the EPIC-PN), and the long 1996 Nov 3 Rosat data. The 1998 Rosat and the Chandra observations are valuable by nature of their several day durations.

TIME SERIES ANALYSIS
Our data originate from instrumentation with widely differing sensitivities: typical count rates vary from 1 count every ∼ 3 s for Rosat HRI to ∼ 6 counts/s for XMM-Newton PN. However, none of these count rates is sufficiently high for a normal distribution of counts to be expected in a time bin of adequate time resolution. The body of literature addressing discrete Fourier Transforms (for example Deeming 1975, Scargle 1982, Schwarzenberg-Czerny 1998 is not directly applicable for our analysis. The treatment of sparse data and event list data is generally accomplished using the Rayleigh Transform (for example de Jager 1991, Mardia 1972. It should be noted that this transform is numerically the same as the discrete Fourier Transform. It will be important in our analysis to define confidence intervals to the derived quantities, in particular the period. In the case where sufficient counts per sample are available, this is generally achieved by least squares fits to the data (for example Kurtz 1985), but this is not possible for sparse or event list data. In this case it is necessary to use the more general maximum likelihood technique, which makes no assumptions on the data distribution (Cash 1979, Lampton, Margon & Bowyer 1976. Using this, Cash (1979) defined the C-and ∆C-statistics. The ∆C-statistic (Cash 1979) has the considerable advantage of being distributed as χ 2 , so that the required confidence intervals can be identified directly from the plot.
The maximum likelihood technique depends on a model probability distribution, which must be specified a priori. Haberl et al. (1997) and Cropper et al. (2001) find the flux variation of RXJ0720.4-3125 to be almost sinusoidal, so in our particular case we can take a sinusoidal variation superimposed on a mean value as an appropriate model. The single frequency and absence of harmonics in the Rayleigh transform of the X00a PN dataset shown in Figure 1 provides the justification explicitly for this approach. The maximum likelihood is then computed simply as equation (35) of Bai (1992). The normalisation of the model is unimportant, as it results in a constant term in the C-statistic, which then is eliminated in the ∆C-statistic. In our case we have two free parameters: normalised amplitude and phase. In general, when performing maximum likelihood fits, the best fit free parameters have to be determined by whatever method may be most appropriate. However, in the case of sinusoids, as noted by Bai (1992), their orthogonality allows them to be determined directly from the Rayleigh Transform (or Lomb-Scargle periodogramme).
The computation of the ∆C-statistic is therefore only slightly more complicated than the periodogramme. It starts in the same way, by a scan through period, calculating the Rayleigh Transform to provide the amplitude and phase at each period. Then at each period the maximum likelihood is calculated using equation (35) of Bai (1992). Finally from equations (4) and (7) of Cash (1979), ∆C is simply twice the difference between the maximum likelihood at a particular period, and the maximum value over the whole frequency range. This results in an inverted "periodogramme", with minimum value of zero. The permitted period range is that between the confidence levels appropriate for χ 2 with the appropriate number of degrees of freedom (1 or 2 depending on whetherṖ is included as a parameter -see section 4).
Details of the procedure are given in Appendix A.
With a period of 8.391 s (Haberl et al. 1997), more than 10 7 cycles have elapsed over the time span of our data set. It is therefore essential to incorporate aṖ term in our analysis techniques when combining datasets at different epochs. The search then becomes a (P0 ,Ṗ ) search, where P0 is the period at the start of the first dataset (the 1993 Rosat observations), and the instantaneous period P = P0 +Ṗ (t−t0) is computed at each point in the grid. Here t0 is the time at the start of the R93 run.
As a more general comment, we recommend this maximum likelihood periodogramme (hereafter MLP) for time series analysis. The uncertainty in a period determination is not evident directly from the width of a peak in a power spectrum computed by classical discrete Fourier transform, Lomb-Scargle periodogramme or Rayleigh transform, as this is set by the window function (c.f. Scargle 1982 appendix D). However, it is directly available from the MLP ∆C statistic: the inverted peak from a sinusoidal signal in data with high signal-to-noise ratio will be narrower than the peak from low signal-to-noise data and the χ 2 can be read directly from the vertical axis (see e.g. Figure 2). Furthermore, as noted above, the MLP technique requires an appropriate model, which may be unknown a priori. In this case the signal can be reconstructed from the harmonic content of the Rayleigh transform or Lomb-Scargle periodogramme, or from the phase folded data, and an appropriate model generated. This model can then be used in the calculation of the MLP, incorporating all of the harmonic information explicitly.

MEASURE OF THE SPIN PERIOD
The only measure of P andṖ presented in the literature for RXJ0720.4-3125 was given by Haberl et al. (1997) using Rosat data. The period was determined to be 8.39115 ± 0.00002 seconds from the 1996 Nov 3 dataset. Based on a linear fit of the data spanning 3 years, these authors derived a 90 % confidence range of −6.0 × 10 −12 to 0.8 × 10 −12 foṙ P , compatible with no period change.
It is computationally unfeasible to explore this range in the P0Ṗ plane in a single analysis of the entire dataset to determine refined values for P0 andṖ . A reduction in the range is required. This can be accomplished initially by extending the Haberl et al. (1997) analysis to include the XMM-Newton data. This is sufficiently accurate to identify the correct peak in the alias patterns in the extended duration R98 and Ch00 data, leading to a further limiting of the permitted range in the PṖ plane. The next stage is a grid search in P andṖ over this restricted range for the 1993 and 1996 combined Rosat data only. The final stage is a grid search using all data in the restricted zone for which there are acceptable solutions from this Rosat analysis.

Individual datasets
We begin by performing an MLP assumingṖ = 0 on each of the longer pointing (R93, R96d, X00a, X00b) in the period range (8.375-8.400) s, incorporating the value of 8.39115 s found by Haberl et al. (1997). The results are shown in Figure 2 (left). As we can see, there is no ambiguity in the period determinations from the R96d, X00a and X00b datasets; the best fit frequencies and their uncertainties are reported in Table 2. A linear least square fit using the 68% formal errors from the MLP results in a value of P0 = 8.39113 ± 0.00011 s, P = 0.0 ± 5.5 × 10 −13 s/s, where, as throughout, P0 is referenced to the start of the R93 run.
The proximity in time between the Ch00 and X00a observations, together with the above upper limit inṖ , permits an unambiguous determination of the second highest peak in the Ch00 power spectrum (Figure 2, right). The same is true for the R98 observations, in which the highest peak is selected. Adding these to the linear least squares fit results in a value of P0 = 8.39107 ± 0.00005 anḋ P = 2.7 × 10 −13 ± 2.5 × 10 −13 . As may be expected, the P0 andṖ are highly correlated. The 68, 90 and 99% confidence intervals are shown in Figure 3. The duration and high count-rate of the X00a PN run alone provides a strong constraint in the (P0 ,Ṗ ) plane. The MLP 68% and 90% confidence intervals are also shown in Figure 3. This indicates that acceptable values for P0 anḋ P lie within the range 8.39106 to 8.39115 s and −1 to +2 × 10 −13 s/s, respectively. We checked the effect of a 10 −7 error in the XMM-Newton clock calibration, as reported recently (Kuster, M., 2001, reported at "New Visions of the X-ray Universe in the XMM-Newton and Chandra Era"). The change in the confidence region is negligible.

Combined Rosat dataset
The P0Ṗ range identified in Figure 3 is sufficiently restricted for an MLP to be performed on the combined 1993 and 1996 Rosat datasets. The use of combined datasets requires coherent phasing to be maintained over the whole dataset (coherent phasing was required only within each of the individual datasets used until this point). The (P0 ,Ṗ ) plane was searched over between 8.39095 ≤ P0 ≤ 8.39120 s, and −4 × 10 −13 ≤Ṗ ≤ 8 × 10 −13 s/s, exceeding the parameter range for the 90% confidence limits from the linear least squares fit. The resulting 68% and 90% confidence contours break up into small regions (aliases) distributed along lines in the (P0 ,Ṗ ) plane. These are also shown in Figure 3 As is evident, there is overlap between the Rosat 90% confidence interval contours and the 99% contours of the linear least squares fit to the individual data (but not quite at the 90% level). This overlap region is also consistent with the X00a PN confidence intervals.

Full dataset
With the further restriction available from the fit to the Rosat 1993 and 1996 data subset in Figure 3, a coherent MLP search was made to all of the data (excluding the uncertain R98 run). The search was made for 8.391050 ≤ P0 ≤ 8.391150 s and −3 × 10 −13 ≤Ṗ ≤ 3 × 10 −13 s/s within the parallel lines P0 = aṖ + b defined by intercepts  (Table 3). All other periods can be excluded at the 95.4% level, but two other possibilites cannot be excluded at the 99% level: these are also given in Table  3 for completeness. The contours for the two best fits are given in Figure 4.
We have folded the data on both P0Ṗ solutions (1) and (2) in Table 3 to check that the relative phasing of all individual runs is correct. We show the folded data on these periods in Figure 5. The data from runs R96a and R96b have been combined to increase the signal-to-noise ratio, and we have not shown the data from R98 because of the uncertainty in the time reference for this run (section 2).  Table 3. The four best fit (P 0 ,Ṗ ) values. Refer to Figure 4 for the confidence interval contours. ∆χ 2 is the difference between the χ 2 of a given solution, and that of solution (1). Fits (3) and (4) are formally excluded at the 95.4% level.
Now performing the a postiori check with the BeppoSax data, we find that these data phase correctly with the main datasets for P0Ṗ (1), but not (2). This is shown in the lowest panel of Figure 5 in each case. This allows us to select P0Ṗ (1) (2) Figure 4. The 68% and 90% MLP contours for P 0 andṖ to the complete dataset except for the R98 data, for the two acceptable (P 0 ,Ṗ ) values (1) and (2). Note that the contours in (2) are referenced to the best global fit derived for (1).
(1) as the most likely timing parameters for RX J0720.4-3125. In any case, both acceptable fits to the data have 3 × 10 −14 <Ṗ < 6 × 10 −14 . For the purposes of our further discussion, the difference between theṖ in (1) and (2) is not important.

DISCUSSION
The refined value ofṖ reported in this paper is consistent with, but two orders of magnitude lower than the extrema of the range reported by Haberl et al. (1997). The first implication is that RXJ0720.4-3125 is unlikely to be spinning down under propeller torques. If the propeller is due to the torque exerted by the interstellar medium, it would be (e.g. Colpi et al., 1998) Pprop ≈ 10 −8 n 9/13 v −27/13 10 B 8/13 12 P 21/13 where n is the external density in cm −3 (n ∼ 1 for the interstellar medium), B12 = B/(10 12 G) and v10 is the star's velocity normalized to 10 km/s. By using the values of P andṖ of RXJ0720.4-3125, we obtain B12 ≈ 16n −9/8 v 27/8 10 . Propeller spin-down dominates over dipolar losses if B12 < 25 √ nv −3/2 10 and this, combined with the previous equation, constrains the star's velocity to extremely small values: v10 < ∼ n 1/3 . Moreover, this scenario makes the matching of thermal and dynamical time scales more difficult. On the other hand, the fact that the propeller contribution is negligeable supports larger values of the star's velocity.
In the "unified" model suggested by Alpar (2001), the observed spin-down is associated with propeller effect by a fossil disk. In this case, under the assumption that the X-ray luminosity of the source is supplied by energy dissipation (frictional torque), upper and lower bounds onΩ can be derived by eq.(2) of Alpar (2001). For RXJ0720.4-3125 the observed luminosity is L ≈ 2×10 31 (d/100 pc) 2 erg cm −2 s −1 , where d100 = d/100 pc and d is the distance (Haberl et al. 1997). This gives: that in turn translates inṖ = (P 2 /2π)Ω While this scenario is still consistent with the spin-down measured for RBS 1223 (Hambarayan et al., 2002), the value ofṖ reported here for RXJ0720.4-3125 is well below this range. Thus, for RXJ0720.4-3125, the associated energy dissipation cannot, alone, account for the source luminosity (unless the source is at ∼ 5 pc, which is unrealistic). That also make less plausible an interpretation of the hardness ratio profile in terms of a spatially variable absorbing matter, co-rotating in the magnetosphere . The observed behaviour is more probably explained by the angle-dependent properties of the emitted radiation.
On the other hand, the slow spin-down rate of RXJ0720.4-3125 is still considerable. The other plausible mechanism which may account for such large and stable val- Figure 5. The datasets folded on P 0Ṗ (1) (left) and (2) (right) in Table 3. ues ofṖ is rotational loss by emission of magnetic radiation. For a dipolar magnetic field, which gives a present value of the magnetic field of B = 2.13 × 10 13 G † ‡ . Despite this value being below those typically quoted for magnetars (B ≈ 10 14 − 10 15 G), it is still extreme and close to the critical value BQ = 4.41×10 13 G at which quantizing effects start to be important in shaping the atmospheric emission. Furthermore, a scenario in which this source is powered by accretion from the interstellar medium must be ruled out: for the present values of P and B the corotating magnetosphere will prevent the incoming material to penetrate below the Alfven radius. The corresponding spin-down age is which, given the numerous uncertainties, is marginally compatible with that inferred by the cooling curves: a few 10 5 yrs for a surface temperature of ∼ 80 eV (see e.g. Kaminker, Haensel & Yakovlev, 2001, Kaminker, Yakovlev & Gnedin, 2001Schaab et al. 1997Schaab et al. , 1999. The discrepancy is less significant if we notice that: a) the cooling curves are strongly dependent on the neutron star mass (Kaminker, Haensel & Yakovlev, 2001), and b) what we are probably observing in X-rays is a region of limited size which is kept hotter than the average star surface, as inferred by the analysis of the pulse-shape . On the other hand, t sd is representative of the true age of the source only in the case in which the magnetic field remained almost constant during the star evolution. The same condition applies for the validity of the cooling curves mentioned above, which do not include the extra input of energy released in the neutron star in case of field decay. It is therefore of fundamental importance to address the field evolution: the strength of B and its variation during the neutron star history determine, in fact, the actual luminosity, the lifetime and even the nature of the energy loss from the star. A related question is where the field is anchored: in the core, penetrating the whole star or confined in the crust. There are three meachanisms which are typically proposed for inducing field-decay: ambipolar diffusion in the solenoidal or irrotational mode and Hall cascade (Goldreich & Reisenegger 1992, Heyl & Kulkarni 1998, Colpi, Geppert & Page 2000. In reality all the three processes co-exist with different timescales, and each of them may dominate the evolution depending on B, L at any given time. In absence of more detailed computations, Heyl & Kulkarni (1998) and Colpi, Geppert & Page (2000) tentatively isolated the three processes and proposed some simple, phenomenological rules to mimic the evolution in the three regimes. We stress that these descriptions are far from being comprehensive of all the effects that can modify substantially the results: they are used in a first approximation and have the main advantage † Note that a twisted magnetosphere may lead to a reduction up to an order of magnitude in the inferred polar value of the magnetic field. The considerations about the age, however, remain unchanged (Thompson, Lyutikov & Kulkarni, 2001). ‡ Here and in the following we specify the discussion to solution (1) of  Table 4. Predicted source age and primordial field for three different mechanisms of decay, simulated as in Colpi et al. (2000).
The present values of P andṖ are those of solution (1) in table 3. In all cases, the source is assumed to be born with P = 1 ms.
of having a simple analytical form. By using the expressions of Colpi et al. (2000), we have estimated the source age and the value of the magnetic field at the birth of the neutron star, B0. Results are shown in Table 4. In all cases the star is assumed to be born with a period of 1 ms: results are not strongly dependent on this exact value, provided it is far less than the present period.
As we can see, allowing for a mechanism involving very fast decay, such as the Hall cascade, we find that the source is now ∼ 4 × 10 4 yr old, and is born with a superstrong field B0 ∼ 10 15 G. Such a young age is only marginally compatible with the absence of a remnant and, more important, is not compatible with the low observed X-ray luminosity (as we compared using not only the standard cooling curves mentioned above, but also some cooling curves kindly provided by Geppert & Colpi, private communication; these latter are computed allowing for the extra-heating due to Bdecay from Hall cascade and predict larger luminosities than the former, making the discrepancy even higher). Underluminous models have been recently presented by Kaminker, Haensel & Yakovlev (2001), who accounted for the enhanced neutrino cooling in presence of strong neutron superfluidity. These solutions may match an age of ∼ 10 4 yrs for RXJ0720.4-3125, but, as discussed by the same authors, they must probably be rejected since they fail in the comparison with observational data of a sample of other neutron stars.
On the other hand, both mechanisms involving ambipolar diffusion predict a magnetic field quite stable over the source lifetime and close to the actual value. Accordingly, the predicted age is ∼ 10 6 years in all cases, close to t sd . Below ∼ 10 14 G the cooling curves are not significantly influenced by decay through ambipolar diffusion (Heil & Kulkarni, 1998), thus, as in the case of constant B discussed above, the scenario is compatible with the observed luminosity. The larger age is also compatible with the absence of a remnant.
If our conclusions are valid, the connection between dim INS and AXPs is not so obvious. RXJ0720.4-3125 has a strong, but not superstrong, field which is compatible with those of the canonical radio-pulsars which have passed the death line. On the other hand, having excluded accretion, what mechanism causes an X-ray emission concentrated in a fraction of ∼ 60% of the star surface remains a mystery, as well as the related question about the validity of using the observed blackbody temperature to locate the source in the cooling diagram. The variation of the surface cooling temperature with the latitude predicted so far for strong fields is not large enough and more complicated explanations are required.

APPENDIX A: DETAILS OF THE MAXIMUM LIKELIHOOD PERIODOGRAMME
The C-statistic is derived from the maximum likelihood ratio and is (Cash 1979), where E is the sum of the expected number of counts according to the model distribution, ni is the observed number of counts in an interval and Ii is the expected number of counts at the time of event ti according to the model distribution. N is the total number of counts. For event list data, the sampling duration tends to zero and ni = 1. The ∆C-statistic is where p is the total number of parameters, and q is the number of free parameters. This is distributed as χ 2 with q degrees of freedom (Cash 1979). (Cmin)p is the minimum C-statistic when all parameters are varied -the minimum over all points in the search grid in the free parameter space. (Cmin) T p−q is the minimum C-statistic at any particular point in the search grid. Substituting, The expected number of counts Ii in equation (A1) is calculated from the model distribution. In our case of a pure sinusoidal variation, because the expected arrival of an event is directly proportional to the model prediction for time ti, we have where A is the normalised amplitude, P is the period and θ0 is the phase. The scaling of this model a0 is unimportant, as constant factors are eliminated when substituting in equation (A3). This model has parameters (A, P, θ0) so p = 3. The periodogramme will scan in P , with best fit values for A and θ0, so q = 1. In the case where we allow a period change, where P0 is the period at ti = 0 andṖ is the period change, the periodogramme will consist of a scan in P0 andṖ . In this case p = 4 and q = 2. In the general case of a maximum likelihood statistic, an optimisation search is required to obtain the best fit for the p − q fitted parameters. Here, in the case of a sinusoid, these can be derived directly from the Rayleigh transform, as pointed out by Bai (1992), where z is the Rayleigh power, and again, ni = 1 for event list data. The amplitude is then This is equivalent numerically to the discrete Fourier transform (DFT) (Deeming (1975), as can be ascertained by reference to Kurtz (1985), equation (A2). The procedure is then for each P (which may be calculated from equation (A5) in a P0 ,Ṗ search) to calculate the amplitude and phase through equations (A7) and (A8) (which amounts to calculating the Rayleigh or discrete Fourier transform), then calculate ∆C through equations (A4) and (A3). As the P grid is searched, the maximum ni ln Ii is stored: ∆C is then simply twice the difference between the calculated value at each period, and this recorded maximum. This difference is distributed as χ 2 for one or two degrees of freedom, depending on whether P is used directly, or calculated in a (P0 ,Ṗ ) search through equation (A5).
An illustration of the power of the MLP is given in Figure 6. Here we have taken the X00a PN and X00b PN data, and computed a discrete Fourier transform in the narrow period range 8.390 to 8.392 s (a much expanded scale by comparison with Figure 2), assuming noṖ . This (upper plot) shows a multiplicity of aliases created by the long gap in the data, superimposed on the broad peak corresponding to the duration of the longest observation (X00a). These are barely resolved in the plot. The lower plot is the MLP with the 68% and 90% confidence levels indicated. This eliminates most of the aliases, as the narrowness of the distribution is related to the precision with which the best fit period can be identified within the broad peak in the upper plot. In the case of the high signal-to-noise ratio X00a and X00b data, this is significantly narrower than the width of the peak, which is set by the window function (Scargle 1982).