Population synthesis studies of isolated neutron stars with magnetic field decay
Abstract
We perform population synthesis studies of different types of neutron stars (thermally emitting isolated neutron stars, normal radio pulsars, magnetars) taking into account the magnetic field decay and using results from the most recent advances in neutron star cooling theory. For the first time, we confront our results with observations using simultaneously the Log N – Log S distribution for nearby isolated neutron stars, the Log N – Log L distribution for magnetars, and the distribution of radio pulsars in the – diagram. For this purpose, we fix a baseline neutron star model (all microphysics input), and other relevant parameters to standard values (velocity distribution, mass spectrum, birth rates …), allowing to vary the initial magnetic field strength. We find that our theoretical model is consistent with all sets of data if the initial magnetic field distribution function follows a lognormal law with and . The typical scenario includes about 10% of neutron stars born as magnetars, significant magnetic field decay during the first million years of a NS life (only about a factor of 2 for low field neutron stars but more than an order of magnitude for magnetars), and a mass distribution function dominated by low mass objects. This model explains satisfactorily all known populations. Evolutionary links between different subclasses may exist, although robust conclusions are not yet possible.
keywords:
stars: neutron — pulsars: general1 Introduction
At the present moment our knowledge about neutron star (NS) evolution is an intriguing puzzle. We know many observational manifestations of young isolated NSs: radio pulsars (PSRs); central compact objects in supernova remnants (CCOs in SNRs); rotating radio transients (RRATs); radioquiet thermally emitting isolated NSs, also known as Xray dim isolated NSs (XDINS) or the Magnificent Seven (M7), and the observational manifestations of magnetars –soft gammaray repeaters (SGRs) and anomalous Xray pulsars (AXPs). Reasons for this apparent diversity as well as possible links between different classes are not entirely clear (see e.g. Manchester et al. (2005); Woods & Thompson (2006); Kaspi (2007); Haberl (2007); Zane (2007); Rea & McLaughlin (2008) for recent reviews about the different subclasses).
In the past few years we have learned that some NSs can show different types of activity, transiting from class to class. For example, PSR J18460258, known for some time as a normal pulsar, demonstrated outbursts typical for AXPs or/and SGRs (Kumar & SafiHarb, 2008; Gavriil et al., 2008). In addition, the total energy release by the object became larger than the rotational energy losses, which according to the original classification discussed in Thompson & Duncan (1996) should place it in the AXP list. Thus, for the first time, we have an example of transformation of a PSR into a magnetar. Several of the SGRs do not show any bursting activity for many years, and if we had not enough information about their violent past, we would have classified them as AXPs. Some of RRATs were shown to emit normal radio pulsar emission (Deneva et al., 2008; McLaughlin et al., 2009). The transient AXP XTE J1810197 and AXP 1E 1547.05408 demonstrated radio pulses (Camilo et al., 2006, 2007). One of the RRATs J18191458 shows thermal properties much similar to the M7 (Reynolds et al., 2006). Hence, divisions between some subpopulations of young isolated NSs (or at least some of their representatives) can be illusive. On the other hand, young, lowfield CCOs (Halpern et al., 2007), normal PSRs ( G, i.e. Crab and Velalike), and SGRs clearly represent NSs born with different properties.
Our brief observational record of NSs ( years at most), low statistics in many cases, and selection effects do not allow to draw a coherent picture of NS “sociology” just from observations. From the theoretical point of view, our understanding of the SN explosion mechanism is not precise enough to provide a solid model of initial parameters of NSs and evolutionary models (thermal, magnetic field, and spin evolution) are related to extremely complicated physical problems (superfluidity and superconductivity in dense matter, electrodynamics in superstrong magnetic fields, etc.) which usually leads to inconclusive results.
In our opinion, it is necessary not only to compare observations with models to verify individual objects, but also to confront theoretical calculations with observational data via population synthesis techniques taking into account as many classes of NSs as possible. Joint constraints by means of simultaneous comparison of theoretical models with different subpopulations should be derived to form a population mosaic. Numerous degrees of freedom in modern evolutionary models must be compensated by different observational tests. Several important studies in this area appeared in recent years, see for example FaucherGiguère & Kaspi (2006); Keane & Kramer (2008), and references therein. In this paper we continue to follow this lane but with an important difference: we attempt to constrain our model and check its consistency using at the same time different populations. This is, to our knowledge, the first time that a multilateral approach is employed. Previous works that focused in a specific NS observational class, ignored that the parameters obtained may be in contradiction with the properties of a different class of objects. For examples, models without magnetic field decay are clearly in contradiction with the existence of magnetars. We try to give a step forward in the direction towards a NS unified model by using at the same time different populations.
Among the different physical ingredients needed to properly model the thermal evolution of NSs, we emphasize that heat transport in the NS crust plays a crucial role during the first million years of its thermal evolution (Yakovlev & Pethick, 2004), which is typically the period during which NSs are detectable with current Xray instruments. Multiwavelength observations in the soft Xray, UV, and optical bands of the thermal emission from a NS’s surface now provide a real opportunity to probe the internal physics of NSs (see Page et al. 2006 for a general overview).
Remarkably, all isolated nearby compact Xray sources that have been detected also in the optical band (RXJ 1856353754, RX J0720.43125, RX J1308.6+2127, and RX J1605.3+3249) have a significant optical excess relative to the extrapolated Xray blackbody emission (Haberl, 2007). An inhomogeneous surface temperature distribution can accommodate this optical excess and can arise naturally if heat conduction in the NS crust is anisotropic due to the presence of a large magnetic field (Geppert et al., 2004, 2006; PérezAzorín et al., 2006). Alternative models are based on anisotropic radiation from magnetized atmospheres, as in Ho et al. (2007), or condensed surfaces as in PérezAzorín et al. (2005). None of the explanations (only temperature anisotropy vs. physical origin) are entirely satisfactory and probably a combination of both effects is needed. Heat conduction can also influence other observable aspects of accreting NSs in low mass Xray binaries, including their quiescent luminosity and the superburst recurrence timescales (Brown et al., 1998). For this reason, one of the goals of this paper is to revisit former population synthesis studies using new NS evolution models that include anisotropic heat transport in NSs crusts and magnetic field evolution (Aguilera et al., 2008a, b; Pons et al., 2009).
In this paper we want to study if the patchy view of
young NSs subpopulations can be explained by a unique set of
smooth distributions of the most important parameters, among which the
the magnetic field distribution plays the main role.
2 Magnetic field decay and cooling model
Very often thermal and magnetorotational evolution of NSs are treated separately. However, in the case of young ( Myr) NSs with magnetic fields G this is incorrect, because temperature affects the electrical resistivity, and therefore the magnetic field evolution, while field decay provides an additional energy source that modifies the temperature of the star. Although for the average radio pulsar population (old and relatively low field) this effect is probably not very important (FaucherGiguère & Kaspi, 2006), there is some observational evidence of the interplay of the magnetic field and temperature during early stages of NS evolution. As discussed in Pons et al. (2007), there is a strong correlation between the inferred magnetic field and the surface temperature in a wide range of magnetic fields: from magnetars ( G), through radioquiet isolated NSs ( G) down to some ordinary PSRs ( G). The main conclusion was that, rather independently from the stellar structure and the matter composition, the correlation can be explained by heating from dissipation of currents in the crust on a timescale of yrs.
This observed correlation has been confirmed later by more detailed 2D cooling simulations combining the insulating effect of strong non radial fields with the additional source of heating due the Ohmic dissipation of the magnetic field in the crustal region (Aguilera et al., 2008a, b; Pons et al., 2009). It was shown that, during the neutrino cooling era and the early stages of the photon cooling era, the feedback between Joule heating and magnetic diffusion is strong, resulting in a faster dissipation of the stronger fields. As a consequence, all neutron stars born with fields over a critical value ( G) reach similar field strengths ( G) at late times. Irrespective of the initial magnetic field strength, the temperature becomes so low after a few million years that the magnetic diffusion timescale becomes longer than the typical ages of PSRs, thus apparently resulting in no dissipation of the magnetic field in old NS. Another interesting result was that the effective temperature of models with strong internal toroidal components is systematically higher than that of models with purely poloidal fields, due to the additional energy reservoir stored in the toroidal field that is gradually released as the field dissipates.
In this paper, we have employed cooling curves obtained from an updated version of the 2D cooling code described in Pons et al. (2009), that also includes the effect of superfluid heat conduction in the inner crust (Aguilera et al., 2009). We have used a Skyrmetype equation of state (EoS) at zero temperature describing both, the NS crust and the liquid core, based on the effective nuclear interaction SLy (Douchin & Haensel, 2001). It is a purely hadronic, relatively stiff EOS that gives typical proper radii in the range km (the radius observed at infinity would be 10–30 % larger depending on the mass), while the crust thickness varies from 0.6 to 1.2 km also depending on the mass of the NS.
We refer to section 4 in Aguilera et al. (2008a) for more details about the cooling models (neutrino emissivities, equation of state, thermal conductivities, etc.). Details about the magnetic field initial geometry can be found in Pons et al. (2009). Our baseline initial model consists of a crustal–confined magnetic field with a poloidal component, parameterized by the value of the radial component at the magnetic pole () combined with a toroidal component with a maximum value of twice (see Eqs. (11) and (13) of Aguilera et al. (2008a)). We find that, although the amplitudes of both fields are of the same order of magnitude, the contribution of the toroidal field to the total magnetic energy is , because this field is non vanishing only in a finite region of the star. This model is in agreement with the results obtained in recent studies of magnetic equilibrium configurations (Ciolfi et al., 2009; Lander & Jones, 2009).
In the left panel of Fig. 1, we show a sample of cooling curves (effective temperature versus true age) for an NS but varying the initial strength of the magnetic field. For G the presence of strong magnetic fields has a visible effect from the very beginning of the evolution (results for a are very similar to this case). The effective temperature of a young, yr magnetar with G is a few times higher than that of a NS with a standard G, and it is kept above K for a much longer time. The effect is qualitatively similar, although somewhat smaller, for high mass stars, and of course it also depends on details about the microphysics input. The influence of the superfluid neutron gap in the core is particularly relevant, in this work we use the results from Baldo et al. (1998). We also use an updated neutron superfluid gap in the crust obtained from Quantum Monte Carlo simulations (Gezerlis & Carlson, 2008). This results in small differences when comparing to the results in Pons et al. (2009). In nonmagnetized NSs, a smaller gap results in higher temperatures at early times (suppression of neutrino emissivity), but varying the gap does not change significantly the temperature of magnetars. The purpose of this paper is to study the observational implications of different initial magnetic fields, being all the rest equal, so hereafter we will fix the underlying physical model and we will only vary the normalization of the field (not its geometry) and the mass of the NS. The variability of the cooling curves on the mass for a fixed magnetic field is shown in the right panel of Fig. 1. We have chosen G and varied the mass of the NS in a wide range, namely and 1.76 . Except for the two lowest masses, the rest of cooling curves are difficult to be distinguished. This is clearly different from the case of nonmagnetized NSs, where there is a clear separation in two scenarios: standard cooling (low mass) and rapid cooling ( in our model) discussed extensively by other authors (Yakovlev & Pethick, 2004; Page et al., 2006). In magnetized NS the fast cooling scenario is masked by magnetic heating, becoming hard to distinguish whether or not a fast neutrino emission process is active in high mass stars, as pointed out in Aguilera et al. (2008a). This raises an important issue: disentangling the magnetic field initial distribution function is needed before we can actually constrain other physical parameters that affect neutron star cooling. For this purpose, we use a population synthesis technique that accounts for joint statistical properties of a set of objects, rather than trying to fit individual objects with particular cooling models.
As in the case of the cooling scenario without extra heating by magnetic field decay, objects with the smallest masses are hotter in average. In the mass distribution function we use they are very abundant and, since they are more easily detectable because of their larger thermal luminosities (for the same field range), most of the observed objects must be low mass NSs. For NSs with the effect of changing mass is controversial, since it strongly depends on the assumptions about neutrino fast cooling processes. Comparison between the left and right panels of Fig.1 shows that changes in magnetic field by a factor of are more important than changing from to for a fixed magnetic field. This implies that population synthesis studies are more sensitive to varying the initial magnetic field distribution than varying the NS mass distribution, unless the population of lowfield ( G) compact objects is discussed. Normally, one would expect that isolated, very cool objects are low field and high mass NSs, so these would be the best candidates to test fast cooling mechanisms without the complications due to the presence of strong magnetic fields.
When we compare our model calculations with observations, we confront emission properties of NSs (observed vs. calculated fluxes or luminosities), not directly their fields, but there is observational evidence for the correlation of these two magnitudes (Pons et al., 2007). In the next section, when we establish constraints on the initial magnetic field distribution, the reader must keep in mind that these are model dependent constraints. What observations actually constrain is the number of luminous objects, which we will translate to magnetic field strength using our theoretical models. As an example, in Fig. 2 we plot the luminosity as a function of the magnetic field strength for a NS model with and at different ages. While the luminosity of young ( yr) objects is less dependent on the magnetic field strength and relatively high (and therefore young NSs are more easily observed), the luminosity of middle aged NSs depends strongly of the magnetic field strength. For a given NS model of a certain age, the luminosity increases sharply above a certain value of the initial magnetic field, thus making more magnetized objects more likely to be observed. We remind again that the axis indicates the initial value of the magnetic field at birth, not the corresponding value at a given age, which is always smaller as discussed above.
3 Log N – Log S distribution for nearby cooling NSs
Previously we performed several calculations of the Log N – Log S distribution for NSs in the solar proximity for different sets of cooling curves. Here we use the model of thermal evolution described above which includes the magnetic field decay. The main motivation is related to the fact that in the standard picture magnetic fields of magnetars decay, and for some of closeby NSs (the M7) there are indications that their fields are G, so additional heating due to field decay can be important. Magnetic fields of some of the M7 sources are estimated by two methods: spectroscopy and spindown rate. The first is based on the unconfirmed hypothesis that wide depressions in their spectra are due to proton cyclotron lines (Haberl et al., 2004). The second method applies the usual magnetodipole braking formula, but this is only possible when a measure of the spin period derivative is available. Estimates due to both methods provide more or less consistent results within a factor of few. In our model with magnetic field decay, the typical strengths of about few G are reached within few yrs for initial fields G, consistent with age estimates for these NSs (see Page et al. (2009) and references therein). Temperatures and spin periods of the M7 sources within our model are also consistent with such ages (see Sec. 6).
To calculate the Log N – Log S distribution (number of objects with a flux above ) of closeby isolated cooling NSs we use the Monte Carlo code developed before (Popov et al., 2003, 2005; Posselt et al., 2008) that builds a synthetic population of nearby isolated NSs. A general description of the population synthesis technique can be found in Popov & Prokhorov (2007). The main ingredients of the present model are: the spatial distribution of NS progenitors, the interstellar medium (ISM) density distribution needed to calculate the observed flux, the NS mass distribution function, a set of cooling curves, discussed in Sec. 2, and the initial magnetic field distribution. We now comment on each of this points before discussing our results.
Progenitors of NSs are distributed according to the model used by Posselt et al. (2008). The contribution due to closeby OBassociations (the Gould Belt) is crucial. As before we consider the region up to 3 kpc from the Sun. Inside 500 pc we take as the distribution of progenitors the distribution of massive stars from HIPPARCOS data (ESA, 1997). Outside this volume, most of NSs (243 per Myr out of the total number of 270 per Myr born inside 3 kpc) originate in one of the OB associations. Others are distributed in the exponential Galactic disc.
In our model for Log N – Log S calculations, the NS formation rate (in units of object per pc per year) is different at different distances from the Sun due to the Gould Belt contribution and the nonuniform distribution of OB associations. Therefore, it is not straightforward to extrapolate the local rate to obtain the total Galactic rate of SNae. Using data from Tammann et al. (1994) we estimate that the value we use ( NS inside 3 kpc per Myr) can be rescaled to 1.2 NSs per 100 years for the whole Galaxy. Here the number 250 is obtained by subtracting the additional contribution due to the Belt ( of all nearby NSs are produced in 600 pc around the Sun). We stress that the value 1.2 per 100 years is only a rough estimate. Sources beyond kpc from the Sun do not contribute much to the Log N – Log S distribution, and in this region the NS formation rate in our model is twice larger than inside 3 kpc.
As we are mostly interested in the Log N – Log S behavior at ROSAT count rates – cts s, the sources at kpc dominate (those born in the associations forming the Gould Belt and in other not very far away OB associations). That is why the global Galactic distribution (Galactic arms, etc.) of NS progenitors is not important for our study with applications to ROSAT cooling NSs.
For the ISM distribution we used an analytical description, which was demonstrated to be successful before (Popov et al., 2003; Posselt et al., 2008). Since here we do not intend to produce an accurate fullsky map of the distribution of sources, and for computational limitations, we do not consider more detailed models for the ISM 3D distribution.
For the mass distribution we use one of the variants, presented in Posselt et al. (2008). The spectrum is derived using HIPPARCOS (ESA, 1997) data about closeby massive stars and calculations by Woosley et al. (2002); Heger et al. (2005). We use eight mass bins (1.1, 1.25, 1.32, 1.4, 1.48, 1.6, 1.7, 1.76 ). The first two bins contribute 30% each. The last two less than 1%. According to this distribution 90% of NSs are born with . Such mass spectrum is in agreement with mass measurements of the secondary components in doubleNS binaries, (see Stairs (2008) and references therein). These NSs never accreted and they can be accepted as good sources for initial mass determination (unless some effects of binary evolution are crucial). In Posselt et al. (2008) it was shown that realistic manipulations with the mass spectrum do not influence Log N – Log S distributions significantly.
We accurately integrate spatial trajectories using the bimodal Maxwellian kick velocity distribution (Arzoumanian et al., 2002) and the potential traditionally used in papers on isolated NSs starting with Paczynski (1990). Nevertheless, the velocity distribution of NSs and the Galactic potential are not important ingredients for the results shown in this section, as we deal with young sources which do not move significantly from their birth places, and their velocities are nearly constant during this time. We tested several velocity distributions including the doubleside exponential with the mean velocity 380 km s proposed by FaucherGiguère & Kaspi (2006). Our results are not sensitive for variation between different velocity distributions which were proposed during last years and which successfully explain the radio pulsar observations. Also, velocity by itself does not influence the observability of sources under study (in contrast with, for example, studies of isolated accreting NSs).
We calculate ROSAT counts for the Log N – Log S plots assuming that the local emission (there is a surface anisotropy consistent with the magnetothermal evolution models) is purely blackbody, i.e. we neglect any nonthermal contribution and we do not consider effects of composition or magnetic fields in the atmosphere. This is partly justified by the fact that the M7 dominate the sample, and their nonthermal emission is negligible (Haberl, 2007). Atmospheric and magnetospheric models can change the spectral energy distribution (not the total luminosity) significantly, resulting in differences in the observed flux at the detectors. However, we think that our sample is too small, and our knowledge about average properties of NS atmospheres is not mature enough to be included in a population synthesis scenario. We suspect that taking this effects into account would not change the results significantly, but in future works, with better knowledge of NSs properties and a larger sample it will be necessary to include these effects explicitly. We must also mention that we considered only heavy element envelopes in our model. The effect of an accreted HHe envelope was discussed in Page (1997) and it will be worth exploring in future works.
Sources are observable in soft Xrays while they are hot. In our previous studies we used cooling tracks till an object cools down to 100 000K. In this paper we put the limit at 300 000K due to computational reasons. However, we tested that this modification does not change our results, as we confront our results with observations of relatively bright sources.
Results of the population synthesis modeling are applied to the ROSAT
allsky survey, which is the most complete survey available at the energy
range of interest (0.11 keV) for our research, and the most uniform sample of closeby
young objects suitable to be used in population synthesis calculations. For
each NS “observed” in the simulation we calculate ROSAT PSPC counts per
second (). Log N – Log S distribution for these sources is shown in
Figs.3, 5.
Filled symbols correspond to additions of one of the M7 NSs,
empty symbols correspond to other sources – Vela, Geminga, B105552,
B0656+14, and 3EG J1835+5918 (the second Geminga).
Note, that for the last five sources we
plot total ROSAT counts (i.e. blackbody plus nonthermal magnetospheric
emission). If we include only the thermal component, then the shape of the
observed Log N – Log S distribution is slightly changed, but clearly the
first and the last filled points will remain as they are, and the change
in the shape is not significant.
Also, nonthermal contributions to luminosities of these PSRs are not large (Becker &
Truemper, 1997).
Error bars represent poissonian
statistical errors (square root of the number of sources). “BSC” is the
upper limit from the ROSAT Bright Source Catalogue (Voges et al., 1999).
In the second Log N – Log S plot we
also show as a horizontal line the most recent limit
at confidence level: for soft, nonvariable NSs,
and for all NSs, including hard and variable sources
(M. Turner et al., in preparation, R. Rutledge, private comm.).
We do not include magnetars in this sample.
If we had considered further regions of the Galaxy
(birth place beyond 3 kpc from the Sun),
we should have included the additional contribution due to young
magnetars, which can be visible from very large distances
(all known young active magnetars are at
distances kpc, according to the McGill SGR/AXP online catalogue
In Fig.3 we show seven curves, each one calculated for a single value of the magnetic field at birth (i.e. all NSs in the modeled population have the same field). This illustrative graph demonstrates that, for our NS model, lowfield NSs ( G) cannot explain the observed sources. If all NSs were born with the same initial field, it should be in the range G. Larger initial fields result in hotter NSs, and therefore in a large number of detectable sources. Here we focus on closeby NSs with fluxes above 0.1 ROSAT PSPC counts per second, as in this range identification of NSs among ROSAT sources is considered to be mostly complete (see, for example, Schwope et al. 1999; Cagnoni et al. 2002 and references therein). This is also confirmed by identification efforts like those by e.g. Agüeros et al. (2006); Chieregato et al. (2005) and references therein. At lower fluxes ( cts s) many sources might be nonidentified, yet. The “second Geminga” (the source with the smallest count rate) is ray selected, and this point must be taken as a lower limit for the Log N – Log S distribution at this flux range. Clearly, multiwavelength studies (cross correlation between catalogues obtained in different bands) are necessary to find more closeby cooling NSs using the existing data, although one of the most crucial issues for identification is probably the Xray positional accuracy. Efforts in this direction, such as in Rutledge et al. (2008), will be very useful.
From results in Fig.3 it is clear that, unlike in previous works, now we have also to worry about the initial magnetic field distribution (Bdistribution), because now our cooling curves depend not only on masses of NSs but also on their fields, and the latter effect seems to be more important. For most of the paper, we have chosen the birth magnetic field to satisfy a lognormal distribution with central value and standard deviation because this type of distribution reproduces the observed distribution of PSRs. The probability of a NS to be born with a magnetic field in the range between and is then
(1) 
where and is the error function.
Model  G  G  G  G  G  G  G  Line  

No mag  0.5  0.5  0.0  0.0  0.0  0.0  0.0  Longdashed  
A1  0.3  0.2  0.1  0.1  0.1  0.1  0.1  Solid  
A2  0.3  0.2  0.2  0.1  0.1  0.1  0.0  Dotted  
G1  1.1  12.5  0.575  0.164  0.114  0.08  0.039  0.019  0.009  Shortdashed 
G2  0.84  13.0  0.37  0.244  0.191  0.126  0.049  0.0165  0.0038  Dotdashed 
G3  0.46  13.5  0.045  0.243  0.396  0.263  0.049  0.0039  0.000075  Dotdotdashed 
In Fig. 4 we show contour plots of the distribution of fits to the Log N – Log S curves in the twoparameter space and , where is given in Gauss. The best fit to the observational data (cross) obtained with the IDL procedure CURVEFIT is and but it is clear that the results of the fit are highly degenerate, and any pair of parameters located in the central diagonal band are allowed. In the figure, we also show the line that corresponds to the family of lognormal distribution functions with an 8% of NSs born with ( G), which interestingly is nearly parallel to the band of lowest . We should stress again that this is a modeldependent result. For our given NS evolution model, this means that actually Log N – Log S curves are not constraining independently the average field at birth or its dispersion (the birth parameters determine the whole evolution), but simply the number of NSs born as highluminosity objects, which according to our underlying physical cooling model is correlated to the number of magnetars. For other cooling models not drastically different from ours, any reasonable field distribution with approximately 810% of NS born as magnetars should also be acceptable. This is consistent with the fact the we do not observe any nearby ( 3 kpc) magnetar. It is important to note that, since magnetars are visible for a long time from large distances, results are so sensitive to the addition of more magnetars that even with 10 observed sources we can place constraints on the fraction of magnetars.
In Fig.5 we present Log N – Log S curves for six Bdistributions (see the Table). Three lognormal distributions (G1, G2, G3), as they are selected from the best fit for Log N – Log S, pass through the observed points. For comparison, we also use several other variants of the Bdistribution, summarized in the Table. Values in the Table correspond to fractions in each magnetic field bin normalized to unity. The first one is an extreme case with no NS with fields above G. The other two are “handmade” distributions (A1 and A2) in both of which 1/2 of NSs belong to the PSRrange fields, and the rest is distributed among highfield objects. The curve for model A2 (dotted line) shows that a 30% of NS with magnetic fields G is already in contradiction with observations of the local population of cooling NSs. Addition of more NSs with very large initial magnetic fields (i.e. model A1) largely overpredicts the number of nearby objects detected as bright thermal sources.
In the presented graphs all Log N – Log S distributions are plotted for 5000 calculated tracks, each of which is used for all eight masses, and all data along a track is used with a time step yrs. So, the results are significantly smoothed. We made a few additional runs for realistic numbers of NSs (810 NSs born during 3 Myrs in the calculated region up to 3 kpc from the Sun). Of course, such Log N – Log S distributions are much more noisy. However, statistical fluctuations cannot change our qualitative results. For example, poissonian error bars for the data points typically bound curves for G1, G2, and G3 distributions. Curves for magnetic field distributions with a significant amount of magnetars cannot explain the data even taking into account statistical fluctuations (unless, of course, a very rare strong fluctuation happened).
Having in mind that observational data can be fitted by different field distributions, and that there are other important parameters not explored (starting from superfluid gaps and ending with properties of the ISM), the important message is that we can reproduce the data on the Log N – Log S plot with realistic distributions, a normal NS model, and without fine tuning. This gives us an opportunity to put a constraint on the number of magnetars. For our magnetothermal evolution model, initial field distributions with more than 30% of NSs with initial fields above G (even if all others have low fields) can be ruled out. Assuming a lognormal initial field distribution, this also implies that the fraction of NSs born with G can hardly be above a few percent. Varying other parameters of the model can modify this conclusion, but not dramatically, unless we allow for significant variations of the birth rate.
4 Log N – Log L distribution for Galactic magnetars
Motivated by a recent important paper by Muno et al. (2008), and as a consistency check of the results of the previous section, we also consider a simple model for the Log N – Log L distribution of highly magnetized NSs. In Muno et al. (2008) the authors analyze a large set of Chandra and XMMNewton data (nearly one thousand exposures) to search for new magnetars looking for pulsating sources in the range 5200 seconds. No new candidate was found, but this fact could be used to place upper limits on the total number of Galactic magnetars with different properties. These authors estimate that the number of magnetars with erg s and pulsed fraction larger than 15% is below 540, and the number of easilydetectable magnetars is . Of course, some (perhaps many) magnetars can have lower pulsed fractions and not being detected as Xray PSRs. In this respect the numbers by Muno et al. (2008) are representative of a fraction of the total population, and total limits maybe larger than claimed. These are, to our knowledge, the best observational limits on the number of such objects. On the other hand, the fact that we observe some magnetars can be used to place lower limits. These limits are too weak to favor one model against another, but it is illustrative to compare some of the models that best fit the closeby, intermediate field NSs to the magnetar population.
In our model the luminosity is of thermal origin, because we do not consider magnetospheric processes that may cause the nonthermal emission. However, since in magnetars both thermal and nonthermal components are supposed to be powered by magnetic field decay, and energy conservation is satisfied by construction in our evolutionary models, we can conclude that the total energy release is calculated correctly, although the spectrum can be different.
We must stress that in this part no Monte Carlo simulation is done. Instead, we use complete cooling tracks of NSs with different masses and magnetic fields to estimate the whole Galactic population of NSs with a given luminosity. Absolute numbers are obtained by normalization to the total birth rate of NSs. We use the Galactic NS formation rate equal to 1/30 yrs. This is close to the upper limit for NS formation rate Keane & Kramer (2008). The uncertainty in the NS formation rate is a factor (see also Keane & Kramer (2008) and references therein) which may shift curves in Fig.6. This NS formation rate is not directly related to the rate used for the Log N – Log S calculations above. In the case of closeby cooling NSs the rate of NS formation is determined by the properties of the Gould Belt and closeby ( kpc) OBassociations. Here, in the case of magnetars, we are interested only in the global galactic rate of NS formation (see a similar discussion in Gill & Heyl 2007).
As we calculate distribution in luminosity, not in flux, we do not take into account interstellar absorption. Calculations are done for the same cooling curves, Bdistributions, and mass distribution as used for Log N – Log S calculations. Formally, very nearby sources with low fields can also contribute to the observed Log N – Log L distribution, but in the range of luminosities we are interested in ( few erg s) their contribution is negligible.
In Fig.6 we show Log N – Log L distributions calculated for the same six Bdistributions described in the Table. We compare our curves with the data by Muno et al. (2008): the upper limit to the number of magnetars with erg s and pulsed fraction larger than 15% is 540 and the number of easilydetectable magnetars is . Since Muno et al. (2008) consider as easily detectable magnetars two types of objects: bright magnetars with small pulsed fraction, and dim, but with very large pulsed fraction, we show this as a rectangular region through which satisfactory models should pass. As for real observations, at the moment we know 5 SGRs (plus candidates), and 10 AXPs (plus transient objects and candidates, see the McGill group online catalogue). These are also shown in the figure (diamonds) for comparison. In the online catalogue the luminosity is given for the range 210 keV. Some magnetars also demonstrate significant hard Xray emission (Mereghetti, 2008). It is not included in the plot, but in logscale the shift is not crucial. Note that this sample is not complete, so it must be taken as a lower limit for the prediction of theoretical models. At the very bright tail it is possible that the sample is close to complete, because there should be no more very bright erg s magnetars.
Modeled Log N – Log L distributions start to flatten at erg s. This value corresponds to the weakest known magnetars. The results from the previous section show that all acceptable lognormal distributions for middleaged thermally emitting NSs are similar in this luminosity range. They all predict about a thousand NSs above this luminosity. Note that we assume that all sources are persistent and the models considered in this work only include steady magnetic field Ohmic decay. In addition to purely Ohmic decay, Pons & Geppert (2007) found that the Hall drift may contribute noticeably to accelerating the dissipation of magnetic fields in young NSs. This Hall phase lasts a few – years and is characterized by an intense exchange of magnetic energy between the poloidal and toroidal components of the field and by the redistribution of magnetic field energy between different scales. It can be expected that such rearrangements and the relatively rapid field decay may enhance the average luminosity and result in crustal breaking and active stages (bursts, flares), as can be observed in magnetars. This would increase the number of bright sources with erg s, and this expectation is preliminary confirmed by some artificial models in which we tried to take into account the possibility of this transient behavior in young highly magnetized NSs. Alternatively, one can also consider a fraction of NSs with larger internal toroidal fields that have larger luminosities, but this introduces yet another free parameter in the problem. A careful study of the first few thousand years of an ultraluminous magnetar’s life and its transient epochs is out of the scope of this paper, that focuses on longterm evolution and statistics.
5 Evolution of PSRs and the diagram
We turn now to PSRs. We have performed Monte Carlo simulations to generate a synthetic PSR population and confront our models with observations. The methodology employed in the simulations closely follows the work by FaucherGiguère & Kaspi (2006) but some parameters of their model are allowed to change according to the results of our previous sections. The main goal of this section is to answer the following question: can we obtain a synthetic PSR population compatible with the observed one and consistent with our previous description for magnetars and closeby isolated neutron stars ? To this end, we start from the optimal population model parameters obtained by FaucherGiguère & Kaspi (2006) and modify only the initial period and magnetic field distributions to account for the effect of magnetic field decay consistent with our model.
To generate the PSR synthetic population we first choose the parameters of
the NS at birth closely following the model described by FaucherGiguère & Kaspi (2006), which we
briefly summarize.
The age of the NS is chosen randomly in the
interval , where 500 Myr. This is shorter than the
age of the Galactic disk, but it is enough for our purpose because PSRs
older than this age have crossed the death line (assuming standard magnetic dipole
braking) and are no longer visible as
PSRs. The place of birth is obtained according to the distribution
of their progenitors (massive Population I stars) which are mainly
populating the Galactic disk and more precisely its arms. The velocity at birth is distributed
according the exponential distribution with a mean value of 380 km s.
The spin period of the star at birth, , is chosen from a normal distribution with a mean value of and standard deviation . Of course, only positive values are allowed. The initial magnetic field at the magnetic pole is obtained from a lognormal distribution with mean value and standard deviation .
Once we have chosen the properties of the NS at birth we solve the appropriate differential equations to obtain the position, period and magnetic field at the present time. We use a smooth model for the Galactic gravitational potential (Kuijken & Gilmore, 1989; Carlberg & Innanen, 1987). The period evolution is obtained by assuming that the rotation energy losses are due to magnetic dipolar emission (orthogonal rotators), where the magnetic field is obtained from our magnetothermal evolutionary models described in section 2.
At the end of the Monte Carlo simulation we end up with a synthetic population of PSRs to be compared with a given observed sample. We use the PSRs detected in the Parkes Multibeam Survey (PMBS) sample (Lyne, 2008) and, to limit the contamination of our sample by recycled PSRs, we further ignore the PSRs with ms or , and those in binary systems. With this restrictions, our resulting sample contains 977 objects. We use the parameters for detectability in the survey, radio luminosity and beaming given in FaucherGiguère & Kaspi (2006). For comparison, The diagram for the sample of real pulsars retained in our analysis is shown in Fig. 7, together with some evolutionary tracks of the model used in the analysis. For lower fields, NSs move nearly along constant magnetic field lines. We must remark again that, in magnetars, a somewhat enhanced field decay is expected to happen during the initial Hall stage. This nonlinear term is not yet included in our numerical simulations and we expect a more vertical initial trajectory for objects in the upper right corner of the diagram.
The problem of defining an optimal model has been well discussed in Sect. 3.7 of FaucherGiguère & Kaspi (2006), and we have adopted their approximate approach that requires of some human judgment, rather than attempting to cover a huge parameter space with our limited computational resources. We have explored the parameter space in the region that best fits the properties of thermally emitting NSs and magnetars, as described in previous sections. Without pretending to perform a rigorous fully quantitative analysis, we have considered the KolmogorovSmirnov goodnessoffit test to quantify our statistical analysis.
In Fig. 8 we show diagrams for typical Monte Carlo realizations and their corresponding distributions of observed PSR periods and magnetic fields (averaged over 50 realizations of 977 detectable PSRs). The upper and lower panels show the optimal models of FaucherGiguère & Kaspi (2006) and this work, respectively. The values for the distribution of the optimal model assuming there is no field decay, are and =0.55, s, and s. Note that FaucherGiguère & Kaspi (2006) use the magnetic field at the equator and the value they give is .
Among all the models analyzed, we find that our optimal model with realistic magnetothermal evolution of NSs corresponds to , =0.6, s, and s. Since our model includes field decay, we find that our average initial magnetic field is about a factor of 2 larger than that of FaucherGiguère & Kaspi (2006), and the distribution is also slightly wider. Because of this larger average field, the initial period distribution has to be shifted to lower values to obtain a visible synthetic PSR population statistically similar to the observed population of PSRs.
Our main conclusion is that, within the parameter space that best fits the observed population of nearby, thermally emitting NSs and magnetars, we can also find an optimal parameterization that satisfactorily explains the observed PSR population. Therefore, it is possible to describe simultaneously different families of the neutron star zoo with a single underlying physical model. Interestingly, a combined statistical analysis of PSRs and thermally emitting NSs allows to break the degeneracy in the parameter space that arises when we try to work with a single family. Models with or , or with wider or narrower distributions (=0.7 or =0.5) give a much worse result in the KolmogorovSmirnov test, obtaining in all these cases Pvalues s. Therefore, having fixed a NS model and initial magnetic field geometry, and within the confidence region obtained in Log N – Log S analysis, we conclude that the observed distribution of radio PSRs is only consistent with values in a narrow vicinity of the optimal model. In the future, as we improve the NS evolutionary models results may change, but the interesting fact is that a combined analysis turns out to be very restrictive and breaks the degeneracy obtained in the study of populations of only nearby thermally emitting NSs or Galactic magnetars. We leave a more extensive study of the influence of these or other parameters (velocity distribution, birth rates, etc.) for future work.
6 Discussion and final remarks.
In this paper we presented a multicomponent population synthesis study. The final goal in this approach would be to make a complete population synthesis with a unique NS physical model that consistently explains all known types of young NSs just varying parameters such as the NS mass, age, or the strength and geometry of the magnetic field. Still, even after more sophisticated theoretical calculations are available, comparison with the data will proceed by steps, confronting each piece of simulated data with observational data of some population, in a similar way to this paper. One reason for engaging multipopulation studies can be illustrated as follows. From Fig. 3 it is visible that the Log N – Log S distribution can be explained by a single field model (with moderate additional heating, in our framework). However, the Log N – Log L distribution for magnetars, as well as the  plot for radio pulsars, of course, cannot be explained by this model for different reasons. In the first case case, because no bright magnetars will be observed; in the latter, because it will be impossible to reproduce the main part of the pulsar population.
To summarize, we believe that the approach we use is well motivated by a necessity to have a natural model without ad hoc assumptions about different subpopulations. Note that here we tested a very “smooth” model, in a sense that we do not put by hand initially distinct populations (magnetars, M7, PSRs, etc.). In each part of our study (Log N – Log S for the M7, Log N – Log L for magnetars, – for PSRs) we model “just NSs” using the same initial magnetic field distribution (which in our model is the main parameter), without specifying unique particular features for a given subpopulation, as it is usually done (see e.g. Popov et al. 2006; Gill & Heyl 2007; Keane & Kramer 2008).
For example, XDINSs are not a separately defined class with initially distinct properties, but they just appear as a population coming out from a smooth unique initial distribution, and with the following features:

The magnetic field in these sources has significantly decayed from larger initial values, so no magnetarlike activity is present;

Due to their large initial fields spin periods are long, so no radio pulsar activity is observed (may be due to narrow beams);

They are still young enough for the magnetic field to be still decaying ( yrs, according to Pons et al. (2009)) and heating the crust, so that they can be observed as relatively bright thermal sources.
In our model a typical M7like source has an initial field of G. The period is sec at a true age of yrs (spindown age yrs). At this time such an object has G and K. The additional energy being input by field decay in these sources can be estimated in each case by multiplying the average magnetic energy density ( erg cm), by the crust volume ( cm) and dividing by the typical Ohmic decay timescale under such conditions (crust density and temperature), which is s. This gives roughly erg s available from magnetic field decay. A fraction of this energy can be radiated in the form of neutrinos, but the energy reservoir to keep the star warm is still very important.
We tested average ages and typical parameters (periods, magnetic fields, etc.) of NSs which contribute to the range of Log N – Log S between 0.1 and 10 cts s, where all observed sources are located. On average modeled NSs with inital fields 10–3 10 G which contribute a lot to this range have ages 2 10 – 5 10 yrs and fields (at the moment of observation) G. Periods are distributed between 220 sec, in correspondence with observed properties of the M7. Of course, some fraction of lower field NSs are also found in this range, in correspondence with observations of cooling radio pulsars (Vela, Geminga etc.). Detailed study of the modeled population in the range 0.110 cts s is in progress, and will be presented elsewhere.
In any population synthesis approach inevitably it is necessary to make simplifications and to neglect some details. Partly, simplifications made in this study are justified by low statistics of known sources, partly by uncertain properties of objects under study. One of the main problems in a population synthesis study is related to possible correlations between different parameters. Some correlations are irrelevant for our purposes (i.e., the correlation between direction of velocity and spin axis) but others can be crucial. For example, in our study we assume that masses, initial magnetic fields and spin periods are not correlated. However, this is very uncertain for all subclasses of NSs. If these three parameters are correlated our results may change. For magnetars it was proposed that they can have massive progenitors (Muno et al., 2006), and normally more massive progenitors are expected to produce massive NSs (Woosley et al., 2002). Thus, we can expect a correlation between field and mass for magnetars but not for other NS classes. Such correlations, valid only for a not well identified part of a population, are very hard to confirm or to rule out. In Fig. 1 one can see that an increase in the value of the magnetic field is much more influential (on the thermal evolution) than changes in mass. That is why we think that moderate massfield correlation in the case of magnetars has very little influence on our results. Other effects like rotation and massloss would complicate the situation even more.
We also neglect all possible effects related to the fact that a significant fraction of NSs are born in binary systems. Potentially, this can be important for magnetars if their large magnetic fields are generated by amplification due to rapid rotation of the protoNS. The precise mechanism of magnetar formation is still unknown. In some models (Popov & Prokhorov, 2006; Bogomazov & Popov, 2009) magnetars are born only in binary systems where a progenitor core was spunup due to accretion or tidal interaction, similar to the main scenario for gammaray burst progenitors.
In our simplified standard scenario (it is clearly visible in Fig.7) the M7 sources do not seem to be descendants of extreme magnetars, but the evolutionary link with some of the AXPs is possible. However, due to numerical limitations and the very different timescales, our model does not include the possibility of enhanced magnetic field dissipation during the fast initial Hall stage, nor transient phenomena that change a quiet Xray emitter into an active bursting source, and back. Understanding the short term violent behavior of young magnetars may help to reconcile even the highest field objects with the M7.
One possibility to test this hypothesis is to look at the velocity distribution of these types of sources. Indeed, velocity is a good invariant on the time scale of 1 Myr. If AXPs are rapidly moving objects, as it was popular to assume some years ago, but M7 have relatively small velocities, then one would conclude that the two populations are not related. Interestingly, there is a recent measure of the transverse velocity of the magnetar XTE J1810197 (Helfand et al., 2007). The measured velocity is slightly below the average for normal young neutron stars, indicating that the mechanism of magnetar birth need not lead to high NS velocities. On the other hand, recently Motch et al. (2009) reported new velocity measurements for the M7 sources, and one of them (RX J1308.6+2127) appears to be a fast object. In the near future we expect that the measurement of proper motions of some of these sources with new observations will contribute to our understanding of their evolutionary link, if any.
In our study we also assumed that the magnetic field structure is similar for all NSs and we only rescale the normalization value for a fixed geometry. This is obviously an arbitrary choice. There are several issues related with the magnetic field geometry that need further investigation. If the magnetic field is not supported by currents located in the crust, but instead by superconducting currents in the core, the field will dissipate on much longer timescales, and the heating mechanism we assume is less important. The amount of energy stored in an internal toroidal field, compared to the corresponding energy of the measured dipolar component is also unknown, and it is an important parameter that can significantly alter the results. We have chosen toroidal fields which maximum amplitude is a factor of 2 the dipole estimate (however, the energy stored on the toroidal field is only about a 10% because it is confined to a small volume), in agreement with recent studies on MHD equilibrium configurations (Ciolfi et al., 2009; Lander & Jones, 2009), but this remains an open question. Some models we tried with toroidal fields one order of magnitude larger produced much hotter objects, and overpredicted the observed number of isolated NSs and magnetars, unless we reduce significantly the birth rate.
In this paper we have not attempted to adjust our cooling curves for low initial magnetic fields in such a way that the local population of normal PSRs with detected thermal emission is perfectly explained. Future more extensive calculations exploring other parameters (superfluid gaps, crust physical properties, etc.) can help to fit better the population of thermally emitting local NSs, and a better knowledge of the mass spectrum can be important for lowfield stars, too. For now, we can work with the known 7 XDINS (for which field decay is probably important) and 45 normal closeby PSRs (i.e. lower field sources without additional heating). In fact, the lowest curve in Fig.5 predicts only 12 objects with flux larger than 0.1 cts s, while we know four nearby PSRs (Vela, Geminga, B1055–52, and B0656+14). The low field cooling curves are more sensitive to details of the equation of state, superfluid gaps, and neutrino emissivities (and the NS mass) than high field models, and changes of this parameters (i.e. the neutron gap in the core) can shift up or down the curves. Ideally, the low field population should be used in a separate analysis to constrain parameters of the interior physics but, given the problem of very low statistics we are facing, it seems meaningless to split our populations in more subgroups at this stage. We made some tests to see if it can be done but, with the low statistics we manage and with the present day uncertainties about details of NS cooling, we do not think that an extensive investigation of microphysical parameters can contribute much to our understanding of NS properties. We certainly need more data before making strong cases in favor of particular cooling models.
Despite many attempts (see Agüeros et al. 2006; Chieregato et al. 2005 and references therein) the number of M7like sources is not increasing significantly. More M7like NSs can be found using deep XMMNewton and Chandra observations. A good example is the new source found by Pires et al. (2009). Another possibility is to have a sample of ray sources selected by Fermi/GLAST, but they should be not M7like, but PSRs with radio beams not pointing towards us (similar to the “second Geminga”). The most promising way to increase the number of M7like sources is related to the future eROSITA instrument aboard Russian satellite SpectrumXGamma. New sources are expected to be dimmer, younger, hotter and further away than the seven ROSAT sources (Posselt et al., 2008), and eROSITA (with nontruncated sensitivity at low energies) will be a perfect tool for them. For magnetars there is some hope to have many more candidates due to the MAXI instrument aboard the International Space Station (Nakagawa et al., 2009). Then, with increased statistics, it will be useful to come back to study a separate Log N – Log S distribution of closeby cooling PSRs.
To summarize, the methodology we employ in this article looks very promising to constrain NS properties as more data coming from future missions arrive and more PSRs are found by radio surveys. A future increase of the observational data will allow to perform much more detailed population synthesis and to establish evolutionary links (if they exist) between different classes of isolated NSs, and to understand better the general evolution of these sources.
Acknowledgments
We thank D.N. Aguilera and R. Turolla for interesting discussions. We thank the referee for critical remarks. This research has been supported by the Spanish MEC grant AYA 200767626C0302, the Research Network Program Compstar funded by the ESF, and the Russian Foundation for Basic research (grants 070200961 and 090200032). SP thanks the Universities of Alicante and Padova for hospitality.
Footnotes
 Note for the arXiv version: After the final (second revision) version of this paper has been submitted, an important eprint by Kaplan and van Kerkwijk appeared (arXiv: 0909.5218). In this article the authors give new observational arguments in favor of the large role of the field decay in the evolution of the Magnificent Seven. Many of our conclusions coincide with those given by Kaplan and van Kerkwijk.
 http://www.physics.mcgill.ca/ pulsar/magnetar/main.html
 From the point of view of velocity and initial spatial distributions these assumptions differ from those used for Log N – Log S calculations. However, as we describe in Sec. 3, Log N – Log S calculations are not very sensitive to the velocity distribution and to large scale spatial distribution. The first is due to relatively small ages of studied sources. The second, because in Log N – Log S at significant ROSAT count rates closeby sources dominate. The same can be said about differences in the models for the Galactic potential in – and Log N – Log S calculations.
References
 Agüeros M. A., Anderson S. F., Margon B., Posselt B., Haberl F., Voges W., Annis J., Schneider D. P., Brinkmann J., 2006, AJ , 131, 1740
 Aguilera D. N., Cirigliano V., Pons J. A., Reddy S., Sharma R., 2009, Physical Review Letters, 102, 091101
 Aguilera D. N., Pons J. A., Miralles J. A., 2008a, A&A , 486, 255
 Aguilera D. N., Pons J. A., Miralles J. A., 2008b, ApJL , 673, L167
 Arzoumanian Z., Chernoff D. F., Cordes J. M., 2002, ApJ , 568, 289
 Baldo M., Elgarøy Ø., Engvik L., HjorthJensen M., Schulze H.J., 1998, Phys. Rev. C , 58, 1921
 Becker W., Truemper J., 1997, A&A , 326, 682
 Bogomazov A. I., Popov S. B., 2009, Astronomy Reports, 53, 325
 Brown E. F., Bildsten L., Rutledge R. E., 1998, ApJL , 504, L95+
 Cagnoni I., Elvis M., Kim D. W., Nicastro F., Celotti A., 2002, ApJ , 579, 148
 Camilo F., Ransom S. M., Halpern J. P., Reynolds J., 2007, ApJL , 666, L93
 Camilo F., Ransom S. M., Halpern J. P., Reynolds J., Helfand D. J., Zimmerman N., Sarkissian J., 2006, Nature , 442, 892
 Carlberg R. G., Innanen K. A., 1987, AJ , 94, 666
 Chieregato M., Campana S., Treves A., Moretti A., Mignani R. P., Tagliaferri G., 2005, A&A , 444, 69
 Ciolfi R., Ferrari V., Gualtieri L., Pons J. A., 2009, MNRAS , pp 821–+
 Deneva J. S., Cordes J. M., McLaughlin M. A., Nice D. J., Lorimer D. R., Crawford F., Bhat N. D. R., Camilo F., et al. 2008, ArXiv eprints
 Douchin F., Haensel P., 2001, A&A , 380, 151
 ESA 1997, VizieR Online Data Catalog, 1239, 0
 FaucherGiguère C.A., Kaspi V. M., 2006, ApJ , 643, 332
 Gavriil F. P., Gonzalez M. E., Gotthelf E. V., Kaspi V. M., Livingstone M. A., Woods P. M., 2008, Science, 319, 1802
 Geppert U., Küker M., Page D., 2004, A&A , 426, 267
 Geppert U., Küker M., Page D., 2006, A&A , 457, 937
 Gezerlis A., Carlson J., 2008, Phys. Rev. C , 77, 032801
 Gill R., Heyl J., 2007, MNRAS , 381, 52
 Haberl F., 2007, Astrophys. & Space Sciences, 308, 181
 Haberl F., Motch C., Zavlin V. E., Reinsch K., Gänsicke B. T., Cropper M., Schwope A. D., Turolla R., Zane S., 2004, A&A , 424, 635
 Halpern J. P., Camilo F., Gotthelf E. V., 2007, ApJ , 668, 1154
 Halpern J. P., Gotthelf E. V., Camilo F., Seward F. D., 2007, ApJ , 665, 1304
 Heger A., Woosley S. E., Spruit H. C., 2005, ApJ , 626, 350
 Helfand D. J., Chatterjee S., Brisken W. F., Camilo F., Reynolds J., van Kerkwijk M. H., Halpern J. P., Ransom S. M., 2007, ApJ , 662, 1198
 Ho W. C. G., Kaplan D. L., Chang P., van Adelsberg M., Potekhin A. Y., 2007, MNRAS , 375, 821
 Kaspi V. M., 2007, ApSS , 308, 1
 Keane E. F., Kramer M., 2008, MNRAS , 391, 2009
 Kuijken K., Gilmore G., 1989, MNRAS , 239, 651
 Kumar H. S., SafiHarb S., 2008, ApJL , 678, L43
 Lander S. K., Jones D. I., 2009, MNRAS , 395, 2162
 Lyne A. G., 2008, in Bassa C., Wang Z., Cumming A., Kaspi V. M., eds, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More Vol. 983 of American Institute of Physics Conference Series, Parkes 20cm Multibeam Pulsar Surveys. pp 561–566
 Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ , 129, 1993
 McLaughlin M. A., Lyne A. G., Keane E. F., Kramer M., Miller J. J., Lorimer D. R., Manchester R. N., Camilo F., Stairs I. H., 2009, ArXiv eprints
 Mereghetti S., 2008, A&ARv , 15, 225
 Mirabal N., Halpern J. P., 2001, ApJL , 547, L137
 Motch C., Pires A. M., Haberl F., Schwope A., Zavlin V. E., 2009, A&A , 497, 423
 Muno M. P., Clark J. S., Crowther P. A., Dougherty S. M., de Grijs R., Law C., McMillan S. L. W., Morris M. R., Negueruela I., Pooley D., Portegies Zwart S., YusefZadeh F., 2006, ApJL , 636, L41
 Muno M. P., Gaensler B. M., Nechita A., Miller J. M., Slane P. O., 2008, ApJ , 680, 639
 Nakagawa Y. E., Yoshida A., Yamaoka K., 2009, in Kawai N., Mihara T., Kohama M., Suzuki M., eds, Astrophysics with AllSky XRay Observations, Proceedings of the RIKEN Symposium, held 1012 June, 2008. RIKEN, and JAXA Suzuki Umetaro Hall, RIKEN Wako, Saitama, Japan., p.120 New Magnetar Frontiers with MAXI Survey. pp 120–+
 Paczynski B., 1990, ApJ , 348, 485
 Page D., 1997, ApJL , 479, L43+
 Page D., Geppert U., Weber F., 2006, Nuclear Physics A, 777, 497
 Page D., Lattimer J. M., Prakash M., Steiner A. W., 2009, ArXiv eprints
 PérezAzorín J. F., Miralles J. A., Pons J. A., 2005, A&A , 433, 275
 PérezAzorín J. F., Miralles J. A., Pons J. A., 2006, A&A , 451, 1009
 Pires A. M., Motch C., Turolla R., Treves A., Popov S. B., 2009, A&A , 498, 233
 Pons J. A., Geppert U., 2007, A&A , 470, 303
 Pons J. A., Link B., Miralles J. A., Geppert U., 2007, Physical Review Letters, 98, 071101
 Pons J. A., Miralles J. A., Geppert U., 2009, A&A , 496, 207
 Popov S. B., Colpi M., Prokhorov M. E., Treves A., Turolla R., 2003, A&A , 406, 111
 Popov S. B., Prokhorov M. E., 2006, MNRAS , 367, 732
 Popov S. B., Prokhorov M. E., 2007, Physics Uspekhi, 50, 1123
 Popov S. B., Turolla R., Possenti A., 2006, MNRAS , 369, L23
 Popov S. B., Turolla R., Prokhorov M. E., Colpi M., Treves A., 2005, ApSS , 299, 117
 Posselt B., Popov S. B., Haberl F., Trümper J., Turolla R., Neuhäuser R., 2008, A&A , 482, 617
 Rea N., McLaughlin M., 2008, in Yuan Y.F., Li X.D., Lai D., eds, Astrophysics of Compact Objects Vol. 968 of American Institute of Physics Conference Series, Rotating RAdio Transients: multiwavelength observations. pp 151–158
 Reynolds S. P., Borkowski K. J., Gaensler B. M., Rea N., McLaughlin M., Possenti A., Israel G., Burgay M., Camilo F., Chatterjee S., Kramer M., Lyne A., Stairs I., 2006, ApJL , 639, L71
 Rutledge R. E., Fox D. B., Shevchuk A. H., 2008, ApJ , 672, 1137
 Schwope A. D., Hasinger G., Schwarz R., Haberl F., Schmidt M., 1999, A&A , 341, L51
 Stairs I. H., 2008, in Bassa C., Wang Z., Cumming A., Kaspi V. M., eds, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More Vol. 983 of American Institute of Physics Conference Series, Double Neutron Star Binaries. pp 424–432
 Tammann G. A., Loeffler W., Schroeder A., 1994, ApJS , 92, 487
 Thompson C., Duncan R. C., 1996, ApJ , 473, 322
 Voges W., Aschenbach B., Boller T., Bräuninger H., Briel U., Burkert W., Dennerl K., Englhauser J., et al. 1999, A&A , 349, 389
 Woods P. M., Thompson C., 2006, Soft gamma repeaters and anomalous Xray pulsars: magnetar candidates. Cambridge University Press, pp 547–586
 Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
 Yakovlev D. G., Pethick C. J., 2004, ARAA , 42, 169
 Zane S., 2007, ApSS , 308, 259